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CHAPTER 1 
INTRODUCTION 


In this chapter the motivation for the current research project will be given. Furthermore, a 
detailed description of the remaining chapters of this report will be presented. 

1.1 Motivation 


NASA Glenn Research Center has an ongoing research program to develop new technologies 
to improve aircraft engine fan containment systems. The program contains a feasibility study to 
replace metallic containment systems with hardwall containment systems composed of polymer 
matrix composites. In such an application, the composite would be loaded at strain rates up to 
several hundred per second. In designing a composite containment system, the ability to 
correctly model the deformation and failure behavior of the material under high strain rate 
loading conditions is of critical importance. 

Experimental techniques to characterize the behavior of polymer matrix composites under 
low strain rate loading conditions have been well established for many years. Furthermore, 
numerous analytical methods have been developed to model the deformation and failure 
behavior of composites under quasi-static loads. However, test methods and analytical 
procedures for characterizing and modeling these materials under high strain rate conditions are 
still under development. 

As will be discussed further in the background section of this report, experimental results to 
date indicate that the deformation and failure behavior of polymer matrix composites is strain 
rate dependent, particularly at higher strain rates. Furthermore, there are strong indications that 
for carbon fiber reinforced composites, the matrix constituent primarily is responsible for the rate 
dependence of the material. In addition, for high strain rate impact applications, ductile polymers 
are more likely to be used in the composite due to their energy absorbing capabilities. As a 
result, the deformation response of the composite as a whole is nonlinear. Analytical methods 
designed to predict the deformation and failure behavior of polymer matrix composites subject to 
impact loads must then have the capability of properly capturing the rate sensitivity and 
nonlinearities which are present in the material response. 

1,2 Overview of Project 

The goal of this project was to develop a nonlinear, strain rate dependent deformation and 
strength model for the analysis of polymer matrix composites. Since the polymer matrix drives 
both the nonlinearity and rate dependence of the composite, constitutive equations using state 
variables were developed to simulate the deformation response of the matrix constituent. The 
constitutive equations were then implemented into a mechanics of materials based 
micromechanics approach to predict the nonlinear, rate dependent deformation response of the 
composite as a whole. A quadratic failure model was utilized to predict ply ultimate strengths. 
Local failure mechanisms were approximated by appropriate combinations of macroscopic ply 
stresses. Finally, the deformation model was implemented in a transient dynamic finite element 
code to model the experimental results. This step provided the ability to analyze composite 
structures subject to impact loads. 
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Chapter 2 of this report gives a relatively complete general background describing previous 
investigations into characterizing and modeling the rate dependent response of polymers and 
polymer matrix composites. Specifically, experimental efforts to determine the effects of strain 
rate on the properties and response of polymer matrix composites will be described. Analytical 
methods that have been developed to predict the nonlinear deformation of polymers, including 
the effects of hydrostatic pressure, will be discussed. Macromechanical and micromechanical 
techniques that have been used to simulate the nonlinear, rate dependent deformation response of 
polymer matrix composites will be presented. Methods that have been developed to predict ply 
ultimate strengths, varying from very simple to more complex, will also be described. The 
correlation between the previous research and the current study will also be presented. 

In Chapter 3, the constitutive equations that were used to model the nonlinear, rate dependent 
deformation response of the polymer matrix are described. First, an overview of the state 
variable modeling method will be given. The one-dimensional form of the constitutive equations 
will be presented, along with procedures for determining the material constants. The extension of 
the equations to three dimensions will be discussed, and the numerical techniques utilized to 
integrate the equations will be presented. Correlation studies conducted using two representative 
polymeric materials, Fiberite 977-2 (a toughened epoxy) and PEEK (a thermoplastic), will be 
discussed. 

The implementation of the polymer constitutive equations into a mechanics of materials 
based composite micromechanics technique will be described in Chapter 4. An overview of the 
method, including assumptions, will be presented. The specific equations used to compute the 
local and ply level stresses based on total strains and fiber and matrix constitutive properties will 
be derived and discussed. Furthermore, equations to compute the effective inelastic strains in the 
composite ply will be given. The numerical techniques utilized to implement the equations will 
be presented. Verification studies conducted using two representative polymer matrix 
composites, IM7/977-2 (IM7 carbon fibers in a 977-2 toughened epoxy matrix) and AS4/PEEK 
(AS4 carbon fibers in a PEEK thermoplastic matrix), will be described. 

In Chapter 5 the equations used to predict ply ultimate strengths will be discussed. An 
overview of the chosen failure criteria will be given. The ply level failure stresses for the same 
representative composites as above will then be predicted for a variety of fiber orientations and 
strain rates and compared to experimental values. 

The implementation of the rate dependent, nonlinear deformation model into the LS-DYNA 
transient dynamic finite element code will be discussed in Chapter 6. An overview of the finite 
element code and the techniques required to implement a user defined material model will be 
presented. The transformation of the polymer constitutive equations into an incremental format, 
required for LS-DYNA, will be described. The implementation of the composite 
micromechanics equations into LS-DYNA will also be discussed. Verification studies conducted 
using the same materials used in earlier chapters will be presented and discussed. 
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CHAPTER 2 
BACKGROUND 


An overview is now given of the work by researchers in the topic areas investigated in this 
study. Specifically, previous work in determining the effects of strain rate on the material 
properties of polymer matrix composites will be described. Next, methods used to model the 
inelastic rate-dependent deformation response of polymers will be considered, including the 
effects of hydrostatic stresses on the yield response. Methods used to simulate the rate-dependent 
deformation of polymer matrix composites will be presented, and techniques used to predict the 
ultimate strength of composites will be discussed. Finally, the relationship between the previous 
work and the current research study will be stated. 

2.1 Strain Rate Effect on Material Properties 

Several experimental studies have been performed with the goal of determining the effects of 
strain rate on the material properties and response of polymer matrix composites. One method of 
performing such studies involved using the split Hopkinson bar technique, which was utilized by 
researchers such as Harding and Welsh [1], Staab and Gilat [2] and Choe, Finch and Vinson [3], 
The technique has been extensively used in characterizing metals at high strain rates [4], The 
basic technique involves propelling a striker bar into a pressure bar. The pressure bar transmits a 
compression stress wave that propagates through the specimen, which is sandwiched between the 
pressure bar and a transmission bar. Gages on the pressure and transmission bars are then used to 
compute the force and displacement in the specimen based on the wave propagation profiles. 

Harding and Welsh [1] created a modified split Hopkinson bar to allow for the tensile testing 
of unidirectional graphite/epoxy composites at high strain rates. Tests were conducted at strain 
rates ranging from 5xl0' 4 /sec to 450 /sec on composites with a [0°] fiber orientation. The stress- 
strain curves obtained were nearly linear until failure, and there was minimal change in the 
elastic modulus and fracture strength with strain rate. Since the response of a [0°] composite is 
primarily fiber dominated, these results indicated that the graphite fibers in tension have minimal 
strain rate dependence. Harding also conducted split Hopkinson bar tests on glass/epoxy and 
glass/polyester woven composites loaded by using a punch [5]. For both types ot materials 
tested, the study found that increasing the punch speed resulted in an increase in the maximum 
punch load present in the specimen, as well as an increase in the shear strength of the specimen. 

In compressive split Hopkinson bar tests that Choe, Finch and Vinson [3] conducted, the 
ultimate stress and modulus were found to increase with strain rate for unidirectional [0°] 
graphite/epoxy specimens, which was different from what was observed by Harding for tensile 
loading [1], The ultimate stress also increased with strain rate for quasi-isotropic graphite/epoxy 
specimens. These results indicated that strain rate dependence may also depend on whether 
tensile or compressive loads are being applied. 

Staab and Gilat [2] conducted tensile split Hopkinson bar tests and static tests on glass/epoxy 
composites at strain rates ranging from lxlO' 4 /sec to approximately 1000 /sec. Laminates with 
fiber orientations ranging from [±15] s to [±75] s were examined. For each of the laminate 
orientations considered, increasing the strain rate increased the initial modulus, maximum stress 
and failure strain. The difference between the static results and the high strain rate results 
increased significantly as the laminate orientation angle decreased. These results differed from 
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what was observed for graphite/epoxy composites, where the properties of laminates with fiber 
dominated orientations displayed minimal strain rate dependence. These results indicated that the 
properties of the glass fibers varied significantly with strain rate. 

Another type of test technique utilized to determine the response of polymer matrix 
composites at high strain rates involved subjecting the composite to explosive pressure pulse 
loading. Such methods were utilized by Daniel, Hsaio and Cordes [6], Daniel, Hamilton and 
LaBedz [7] and Al-Salehi, Al-Hassani and Hinton [8], Daniel, et. al. [6,7] utilized an expanding 
ring, where a thin composite ring was subjected to an explosive pulse loading. In combination 
with static tests, the effects of strain rate on the longitudinal, transverse and shear properties of a 
unidirectional graphite/epoxy composite system were determined. The experiments were 
conducted at strain rates ranging from lxlO' 4 /sec to 500 /sec. In the longitudinal direction, the 
modulus increased slightly, but the strength and failure strain showed almost no variation with 
increasing strain rate. In the transverse direction, on the other hand, the modulus and strength 
increased significantly with increasing strain rate. When the shear properties were examined, the 
modulus and strength again increased noticeably with increasing strain rate. The results in the 
longitudinal direction were similar to what was found by Harding and Welsh [1], with the 
properties showing little variation with strain rate for a carbon fiber reinforced composite loaded 
along the fiber direction. However, since the transverse and shear properties are matrix 
dominated, the fact that these values varied with strain rate indicated that for a graphite/epoxy 
composite, the behavior of the epoxy matrix does vary with strain rate, and drives the rate 
dependence of the composite. 

Using similar explosive pressure pulse testing techniques, Al-Salehi, et. al. [8] found that for 
a filament wound glass-epoxy tube, the burst strength and failure strain increased with increasing 
strain rate, while the elastic modulus displayed minimal variation with strain rate. These results, 
differing from what was observed for the unidirectional graphite/epoxy composites, indicated 
once again that material type and fiber layup can both have significant effects on the rate 
dependence of a composite. 

The correlation between the strain rate dependence of the composite and the rate dependence 
of the matrix observed in graphite/epoxy composites was further explored by Groves, et. al. [9]. 
Tensile and compressive stress-strain curves were generated for an epoxy matrix over a variety 
of strain rates. For both tensile and compressive loads, the modulus was found to increase with 
strain rate, with particularly significant increases occurring at strain rates above 10 /sec. Both 
tensile and compressive strengths were also found to increase with strain rate. The results from 
these tests confirmed that strain rate has a significant effect on the properties of the matrix in a 
graphite/epoxy system. 

The overall result from these experiments is that the material properties and response of 
polymer matrix composites were found to vary with strain rate. The fiber material and fiber 
layup affected the nature of the strain rate dependence. For graphite/epoxy composites under 
tensile loading, the strain rate dependence appeared to be primarily driven by the strain rate 
dependence of the polymer matrix. These results indicated that in modeling graphite fiber 
reinforced polymer matrix composites, simulating the rate dependence of the polymer matrix 
correctly is of critical importance. 
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2.2 Polymer Constitutive Modeling 


Polymers are known to have a rate dependent deformation response. Traditionally, for very 
small strain analyses, linear viscoelasticity has been used to simulate the material behavior 
phenomelogically [10]. In linear viscoelastic models, combinations of springs and dashpots have 
been used to capture the rate dependent behavior. For cases where the strains are large enough 
that the response is no longer linear, nonlinear viscoelastic models have been developed. For 
example, in a model developed by Cessna and Sternstein [11], nonlinear dashpots were 
incorporated into the constitutive equations. The rate dependence observed in polymer 
deformation has also been modeled empirically by scaling the yield stress as a function of strain 
rate [12]. 

A more sophisticated technique in polymer constitutive modeling has taken a molecular 
approach. In this method [13], the polymer deformation was assumed to be due to the motion of 
molecular chains over potential energy barriers. The molecular flow was due to applied stress, 
and the internal viscosity was assumed to decrease with increasing stress. The yield stress (the 
point where permanent deformation begins) was defined as the point where the internal viscosity 
decreased to the point where the applied strain rate is equal to the plastic strain rate. Internal 
stresses were also defined [13,14]. These stresses represented the resistance to molecular flow 
that tends to drive the material back towards its original configuration. Another approach to 
polymer deformation assumed that the deformation was due to the unwinding of molecular kinks 
[15]. In both approaches, constitutive equations were developed [13-17]. In these equations, the 
polymer deformation was considered to be a function of parameters such as the activation 
energy, activation volume, molecular radius, molecular angle of rotation, and thermal constants. 
Furthermore, the deformation was assumed to be a function of state variables that represented the 
resistance to molecular flow caused by a variety of mechanisms. The state variable values 
evolved with stress, inelastic strain and inelastic strain rate. 

An alternative approach to the constitutive modeling of polymers has utilized, either directly 
or with some modifications, viscoplastic constitutive equations which have been developed for 
metals. For example, Bordonaro [18] modified the Viscoplasticity Theory Based on Overstress 
developed by Krempl [19]. In Bordonaro’s model, the original theory was modified to attempt to 
account for phenomena encountered in polymer deformation that are not present in metals. For 
example, polymers behave differently from metals under conditions such as creep, relaxation and 
unloading. Other authors, such as Valisetty and Teply [20] and Zhang and Moore [21], also 
utilized viscoplastic constitutive equations developed to model the deformation of metals to 
analyze polymers. However, in these studies, only uniaxial tensile behavior was analyzed, and no 
attempt was made to consider phenomena such as unloading, creep or relaxation. 

The conclusion to be drawn from the work discussed above is that it is possible to analyze 
the rate dependent response of polymers by simulating the physical deformation mechanisms. 
Furthermore, the physical deformation mechanisms can be modeled by the use of state variables. 
This approach is very similar to what has been used in viscoplastic constitutive equations in the 
analysis of metals. The work completed to date indicates that modeling techniques developed for 
metals can be adapted for use with polymers. However, appropriate modifications must be made 
to the equations and consideration must be given to the range of applicability of the model. 
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2.3 Pressure Dependence of Polymer Ineiastic Deformation 


The hydrostatic stress state has been found to have a significant effect on the yield behavior 
of solid polymers [12,13,22,23]. Increasing the hydrostatic pressure has been found to increase 
the yield stress [13,22,23], The effects of hydrostatic pressure also result in differences in the 
tensile and compressive yield stresses of a polymer. This effect has been incorporated into the 
Maximum Shear Stress (Tresca) and Maximum Distortion Energy (von Mises) yield theories to 
develop yield criteria that have been applied to polymers [12]. In the Maximum Shear Stress 
criteria, the shear yield stress is a linear function of the hydrostatic stress. In the Maximum 
Distortion Energy criteria, the octahedral shear stress at yield is a linear function of the mean 
stress. Ward [13,22,23] also incorporated these concepts into modifying the Eyring based yield 
criteria. Specifically, the hydrostatic pressure was included as an additional term in the equations 
relating the octahedral shear strain at yield to the octahedral shear stress at yield. 

If the effects of hydrostatic pressure on the inelastic deformation of polymers are to be 
incorporated into a state variable constitutive model, appropriate effective stress definitions need 
to be developed. The effective stress is usually defined as being related to the second invariant of 
the deviatoric stress tensor [24], Bordonaro [18] determined that the usual definition did not fit 
the multiaxial tensile and shear data obtained in their study for Nylon and PEEK. As a result, 
several alternative effective stress definitions were developed. One such modification involved 
utilizing total stresses instead of stress deviators in the effective stress terms. Another 
modification consisted of adding a multiple of the mean stress to the product of stress deviators. 
However, neither of these modifications yielded acceptable results when compared to the 
experimental data. These results indicated that the definition of effective stress might need to be 
modified for the analysis of polymers. 

2.4 Rate Dependent Composite Constitutive Modeling 

A variety of methods have been applied to model the rate dependent response of polymer 
matrix composites. As will be described below, models have been developed at both the 
macromechanical (ply) level and the micromechanical (constituent) level. 

2.4. 1 Macromechanical Approaches 

In the macromechanical approach, the composite material is modeled as an anisotropic, 
homogeneous material, without any attention being paid to the individual constituents. For 
example, Weeks and Sun [25] developed a macromechanical, rate dependent constitutive model 
to analyze the nonlinear rate dependent response of thick composite laminates over a variety of 
strain rates. Building upon previous work conducted by Yoon and Sun [26] and Gates and Sun 
[27], the inelastic behavior of a carbon fiber reinforced thermoplastic was simulated through the 
use of a quadratic plastic potential function. A scaling function was defined to model the 
variation of the response due to varying fiber orientation of a single ply. The finite element 
method was utilized to analyze a composite laminate, where a layer of elements was used to 
simulate a single ply. The rate dependence of the deformation response was captured by varying 
the material properties as a function of strain rate. Thiruppukuzhi and Sun [28] later modified 
this technique in order to directly incorporate the rate dependence of the material response into 
the constitutive model. Espinosa, Lu, Dwivedi and Azvattieri [29] utilized a similar type of 
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approach. However, the finite element procedures were specifically reformulated to account for 
the dynamic, large deformation response often seen in high strain rate impact problems. 

Other approaches have been used to model the high strain rate response of polymer matrix 
composites on the macromechanical scale. For instance, O’Donoghue, et. al. [30] assumed a 
linear elastic deformation response, but reformulated the constitutive equations to separate out 
the hydrostatic and deviatoric stresses. The stress separation facilitated implementation of the 
technique into a transient dynamic finite element code. A similar approach developed by Tay, 
Ang and Shim [31] utilized an empirical model that scaled the “dynamic” stresses (i.e. the 
stresses due to dynamic loading) as a function of strain rate. 

2.4.2 Micromechanical Approaches 

Research has also been conducted in simulating the high strain rate deformation response of 
polymer matrix composites through micromechanics approaches. In micromechanics, the 
effective properties and response of the composite are computed based on the properties and 
response of the individual constituents. Several different types of methodologies have been used 
in micromechanics analyses. These techniques have been thoroughly reviewed and discussed in 
works such as [32-36], In general, three types of methodologies have been used. All of these 
approaches were based on analyzing the behavior of a unit cell of the composite. The unit cell is 
the smallest portion of the composite for which the behavior of the unit cell is considered to be 
representative of the response of the composite as a whole. 

The simplest types of micromechanics techniques developed have been mechanics of 
materials based methods, in which various uniform stress and uniform strain assumptions were 
utilized within the composite unit cell to compute the effective properties and response of the 
material. Examples of this type of approach include the traditional “rule of mixtures” equations 
[33], and the simplified micromechanics equations developed by Murthy and Chamis [37], While 
this approach involved a great deal of approximation and simplification, the resulting equations 
were very simple in form, were very easy to implement within a computer code, and were very 
computationally efficient. 

A more sophisticated method to compute the effective properties of composite materials 
involved using continuum mechanics techniques. In this type of approach, the equations of 
continuum mechanics were solved in an average sense within the unit cell. Examples of this 
methodology include the Concentric Cylinders Model [35], the Self Consistent Method [35], the 
Mori-Tanaka Method [38], and the Method of Cells [32], Continuum mechanics methods more 
completely satisfied the field equations of mechanics, resulting in a more accurate representation 
of the physics of the problem, in comparison to mechanics of materials techniques. However, 
these approaches still often resulted in closed form solutions, which permitted reasonable 
implementation and execution of these techniques within a computer code. 

The most accurate micromechanics techniques have been the numerically based methods. In 
this approach, the fiber and matrix were explicitly modeled using either finite elements or 
boundary elements. The effective response of the unit cell was then computed by conducting a 
finite element or boundary element analysis. Examples of this approach can be found in [39,40], 
A numerical technique was also developed by Walker, et. al. [41,42], in which integral equations 
were developed using Fourier series and Green's function approaches. The integral equations 
were then solved using numerical methods. This type of analysis yielded the greatest accuracy. 
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but the execution times required to conduct the analysis on a computer were often quite 
substantial. 

In simulating the high strain rate deformation response of polymer matrix composites 
through micromechanics approaches, finite element and continuum mechanics methods have 
been used. For example, Espinosa, Emore and Xu [43] utilized a finite element approach in 
which the fibers and matrix were explicitly modeled in the finite element mesh. A molecular, 
state variable based polymer constitutive equation, similar to what was described earlier in this 
report, was used to model the inelastic, rate dependent deformation response of the polymer. 
Clements, et. al. [44] applied the Method of Cells [6] to account for the stress wave propagation 
at the constituent level seen in an impact problem. Aidun and Addessio [45] also used the 
Method of Cells to simulate a high strain rate impact problem. For their analysis, they developed 
a nonlinear elastic constitutive model to simulate the response of the polymer matrix. In addition, 
they reformulated the micromechanics equations to separate out the hydrostatic and deviatoric 
stress components to facilitate implementation of the technique into a transient dynamic finite 
element code. 

2.4.3 Summary 

The overall conclusion from this portion of the review is that the high strain rate deformation 
response of composite materials has been modeled using a variety of methods. Furthermore, the 
nonlinear response of the composite was accounted for within the constitutive equations. In 
macromechanical methods, the nonlinearity and rate dependence of the deformation response 
were accounted for at the ply level. In micromechanical approaches, the rate dependence and 
nonlinearity of the polymer matrix was modeled at the constituent level. The homogenization 
techniques then computed the effective deformation response of the composite based on the 
response of the individual constituents. 

2.5 Ply Level Composite Failure Models 

As indicated in reviews such as those conducted by Nahas [46], many failure criteria have 
been developed to predict the ply ultimate strength in polymer matrix composites. No one 
criterion has been found to be superior for all materials and loading conditions. Several “classic” 
criteria have been used, as detailed in texts such as those by Gibson [34] and Herakovich [47] 
and in review papers such as those by Reddy and Pandey [48]. Simple models such as the 
Maximum Stress and Maximum Strain criteria simply compared the stresses (or strains) in each 
coordinate direction to the ultimate values. The drawback to these models was that stress 
interaction was not accounted for to any significant extent. In the Tsai-Hill criteria, on the other 
hand, a quadratic combination of stresses and strengths in each coordinate direction was 
compared to a failure value. In this model, stress interaction was accounted for, but differences in 
tensile and compressive strengths were not considered. The Tsai-Wu equation utilized a tensor 
based failure criterion, which allowed for differences between tensile and compressive failure 
strengths. In general, the more sophisticated methods such as Tsai-Hill or Tsai-Wu provided 
improved failure predictions than simple models such as Maximum Stress or Maximum Strain. 
However, studies such as that conducted by Sun and Quinn [49] showed that this is not always 
the case. 
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In the models discussed above, composite failure was considered on purely a 
macromechanical scale, with no attention being paid to the specific mechanisms which cause 
failure. In reality, failure in a composite is a result of specific local mechanisms such as tensile 
fiber failure, matrix cracking or delamination [33]. While accurate modeling of these local 
failure mechanisms would require a detailed micromechanical analysis, several researchers have 
developed phenomenological failure criteria that predict local failure mechanisms based on ply 
level stresses. 

Hashin [50] developed a failure model that predicts fiber failure or matrix failure (tensile or 
compressive) based on appropriate quadratic combinations of stresses and failure strengths. 
Chang and co-workers [51-53] developed similar quadratic failure criteria, but also accounted for 
nonlinearities in the stress-strain curve by using integrals of strain energy in the shear terms 
instead of shear stresses and shear strengths. While the Chang criteria assumed a plane stress 
condition. Yen [54] extended the Chang criteria to account for a three-dimensional stress state. 
Yen then utilized the modified criteria to predict composite damage and failure for a transient 
dynamic impact application. Banerjee [55] also successfully utilized the Hashin and Chang 
failure criteria to predict damage and failure for a polymer matrix composite subject to impact 
loads. Rotem [56] developed a quadratic failure criteria based on local failure mechanisms, but in 
the matrix failure equations the matrix failure strength and matrix stress were included along 
with the ply level stress and strength values. Tabiei et al. [57] utilized the Hashin criterion to 
predict the failure of a polymer matrix composite with a nonlinear deformation response. 
However, the nonlinear behavior was accounted for in the composite constitutive law instead of 
in the failure criteria. Langlie and Cheng [58] developed failure criteria based on local failure 
mechanisms to be used for composites subject to high strain rate impact. Their criteria used a 
maximum stress type of approach instead of a quadratic equation including stress interaction to 
predict the local failure modes. Pecknold and Rahman [59] applied a micromechanics model to 
predict the deformation response of the composite. The constituent level stresses computed using 
the micromechanics were then compared to the constituent strengths in order to predict failure. 
The research studies discussed here indicate that a failure criterion based on approximating local 
failure mechanisms might produce good predictions. A mechanism based failure criterion would 
also facilitate the development of a material degradation model when applied to the analysis of 
structures composed of composite materials. 

2.6 Relationship of Literature Review to Current Project 

The overall goal of the research study presented in this report was to conduct a preliminary 
investigation into some of the issues involved in modeling high strain rate impact of polymer 
matrix composites. Therefore, the goal was not to develop a refined, sophisticated set of 
analytical methods that could completely simulate an impact event. Instead, the goal was to 
investigate some of the relevant issues involved in the problem, and to see how the different 
pieces of the analysis would fit together to completely solve the problem. As a result, in this 
work, simplified methods, or modifications to previously existing methods, were used to conduct 
the analyses presented here. Once these preliminary investigations have been completed and the 
relevant issues identified, each of the aspects of the analytical model could then be refined and 
adjusted, eventually leading to the development of a complete deformation and failure model for 
simulating the impact of polymer matrix composites. 
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The first step in this process was the development of constitutive equations to model the rate- 
dependent, inelastic deformation response of ductile polymers, as the rate dependence of 
polymer matrix composites has been found to be driven by the matrix. As mentioned in Section 
2.2, state variable methods have been successfully applied to the rate-dependant analysis of 
polymers. For this work, the Ramaswamy-Stouffer state variable constitutive equations, which 
have been used to model the viscoplastic response of metals, was modified to analyze the 
inelastic response of the polymers considered here. 

Section 2.3 described a variety of methods that have been used to simulate the deformation 
response of polymer matrix composites. For this work, since the deformation response of the 
polymer matrix was to be explicitly modeled, micromechanics techniques were used to predict 
the effective composite response. As discussed earlier, a variety of techniques, from fairly simple 
to quite complex, have been developed to conduct micromechanical analyses. For this study, 
since the deformation model was implemented within a finite element code, the desire for this 
preliminary investigation was to keep the analysis method fairly simple while still capturing 
most of the physics of the problem. Utilization of a fairly simple micromechanics approach 
facilitated the finite element implementation. As will be discussed in Chapter 4, a previously 
developed mechanics of materials approach was modified to carry out the micromechanics 
calculations. As discussed in Chapter 6, the micromechanics equations, with the polymer 
constitutive equations embedded within them, were implemented into a commercially available, 
transient dynamic finite element code and the tensile response of representative composites was 
simulated. 

In predicting the failure of polymer matrix composites subject to impact loads, an important 
first step is the prediction of the failure of the individual plies. However, in order to incorporate 
property degradation models into the finite element analysis, the particular mode of failure (fiber 
or matrix failure) must be identified. As will be discussed in Chapter 5, for this preliminary study 
the goal was to predict ply failure based on local failure mechanisms using ply level stresses and 
strengths. The Hashin failure criteria discussed in Section 2.5 appeared to meet the desired 
criteria and was used in this study. 
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CHAPTER 3 

POLYMER CONSTITUTIVE EQUATIONS 


In this chapter, a set of constitutive equations based on the state variable method are 
presented which allow for the analysis of the nonlinear, rate-dependent deformation response of 
ductile polymers. An overview of the state variable modeling method is discussed first. Next, the 
one-dimensional simplification of the constitutive equations is presented, along with the 
procedures for determining the material constants. The equations are then extended to three 
dimensions. The algorithm used to numerically implement the constitutive model is discussed, 
and the equations are characterized and correlated for two representative polymers: Fiberite 
977-2 and PEEK. 

3.1 State Variable Modeling Overview 

As discussed previously, the rate dependence of a polymer matrix is primarily a function of 
the rate dependence of the matrix constituent. Furthermore, the stress-strain response of 
polymers is nonlinear above one or two percent strain [10]. Consequently, a need exists for 
constitutive equations that capture the nonlinear, rate dependent deformation response of the 
matrix material. 

Constitutive equations for polymers that incorporate the deformation mechanisms of the 
material have been developed. For polymers, deformation is due to the motion of the molecular 
chains [13]. At small deformation levels, prior to yield, there is also a resistance to the molecular 
flow. In constitutive equations, a state variable approach has been utilized to model the 
mechanisms that cause material deformation [24]. In this methodology, the specific changes in 
the local details of the material microstructure were not considered. Alternatively, variables were 
defined which were meant to represent the average effects of the deformation mechanisms that 
were present. These variables evolved as a function of external parameters such as the stress, 
inelastic strain, and the current value of the state variable. Furthermore, the inelastic strain rate 
was defined to be a function of the state variables and external variables such as the current 
stress level. 

The state variable approach to constitutive modeling has been fairly extensively utilized to 
model the inelastic deformation response of metals, which exhibit a viscoplastic response above 
about one-half of the melting temperature. An important characteristic of this method was that 
there was no defined yield stress or onset of inelasticity [24], Alternatively, inelastic strain was 
assumed to be present at all values of stress. The inelastic strain was just assumed to be very 
small compared to the elastic strain at low stress levels. When the inelastic strain evolved to a 
more significant value, the stress-strain curve began to exhibit a nonlinear response. 
Furthermore, a single, unified variable was utilized to represent all inelastic strains. The effects 
of viscoelasticity, plasticity and creep were not separated in this type of approach, but were 
combined into one unified strain variable. 

There is some physical motivation to utilize state variable models that were developed for 
metals to simulate the nonlinear deformation response of polymers. For polymers, while the 
nonlinear deformation response is due to the nonlinear response of long chain molecules as 
opposed to the propagation of dislocations for metals, a unified inelastic strain variable can still 
be utilized to simulate the nonlinear behavior. In addition, the “saturation stress” in metals and 
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the “yield stress” in polymers (the point where the stress-strain curve becomes flat) have both 
been defined as the stress level at which the inelastic strain rate equals the applied strain rate in 
constant strain rate tensile tests [13, 24], 

In metals, the inelastic strain rate has often been modeled as being proportional to the 
difference between the deviatoric stress and “back stress” tensors. The back stress was defined as 
a resistance to slip resulting from the interaction of dislocations under a shear stress with a 
barrier. As dislocations pile up at a barrier, atomic forces will cause additional dislocations 
approaching the barrier to be repelled. This repelling force was referred to as the back stress. The 
back stress is in the direction opposite to the local shear stress in uniaxial loading (and thus will 
be orientation dependent in three-dimensional loading). The net stress producing slip or inelastic 
strain is related to the difference between the shear stress and the back stress [24], An isotropic 
initial resistance to slip is also present in metals due to the presence of obstacles such as 
precipitates, grains and point defects. This initial resistance to slip has often been referred to as a 
“drag stress”. Similar concepts have been used in the deformation modeling of polymers. Ward 
[13] discussed how creep strain and plastic strain could be defined as being proportional to the 
difference between the applied stress and an “internal stress”. The internal stress was defined as 
evolving with increasing strain. 

It is important to note several significant limitations in using equations developed for metals 
to simulate the nonlinear deformation response of polymers. Polymers exhibit nonlinear strain 
recovery on unloading, while metals display linear elastic strain recovery. The unloading 
behavior of polymers may not be represented accurately with a constitutive equation that was 
developed for metals. Furthermore, phenomena such as creep, relaxation and high cycle fatigue 
may not be simulated correctly in polymers with a metals based constitutive model. However, for 
predicting failure in high velocity impact loading, the ultimate future goal of this research, none 
of these phenomena were considered to be extremely significant. Furthermore, in the analyses 
described in this report, only uniaxial tensile loading was considered. 

The constitutive equations utilized here will most likely only be valid for relatively ductile 
polymers with at least some level of crystallinity, such as thermoplastics or toughened epoxies 
(with thermoplastic tougheners), where the correlation between the polymer deformation 
response and the deformation response of metals is strongest. However, it is not clear at this time 
exactly how the level of crystallinity in the polymer will affect the applicability of the equations. 
To fully simulate the complete range of polymer deformation response, a combined viscoelastic- 
viscoplastic model would most likely be required. A combined viscoelastic-viscoplastic model 
would have the ability to completely represent all of the mechanisms present in polymer 
deformation. Efforts have been made by several researchers to develop such a model [16,17,60]. 
However, such a combined model, which could be very complex, was not considered for 
this study. 

3.2 One-Dimensional Constitutive Equation 

A state variable based constitutive model was utilized to simulate the rate dependent, 
nonlinear deformation behavior of the polymer matrix. For the purposes of determining material 
constants and model correlation, one-dimensional uniaxial versions of the equations were 
considered first. An important point to note is that the flow equation, in its three-dimensional 
form, was based on deviatoric stresses and stress invariants, which are the primary drivers for the 
inelastic deformation of both polymers and metals [12,24]. Only the equations for the inelastic 
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strain rate are presented. The total strain rate is the sum of the elastic and inelastic strain rates. 
The total strain rate for uniaxial loading was defined as follows: 


e = 



(3.1) 


where £ was the total strain rate, a was the stress rate, e' was the inelastic strain rate, and E 
was the elastic modulus of the material. 

Several key assumptions were made in the constitutive model. First, even though in high 
strain rate impact situations adiabatic heating may be a significant issue, for this study 
temperature effects have been neglected, and all of the results were obtained for room 
temperature. Second, small strain conditions have been assumed. In reality polymers, particularly 
in compression, can be subject to very large strains. Furthermore, structures undergoing high 
strain rate impact are subject to large deformations and rotations. However, incorporating large 
deformation and rotation effects into the constitutive equations would add a level of complexity 
that is beyond the scope of this study. For the purposes of this modeling effort, therefore, all 
large deformation effects have been neglected. Future efforts will include modifying the 
constitutive equations to account for large deformation effects. 

For this study, constitutive equations developed by Ramaswamy and Stouffer [24] were used 
to compute the inelastic strain rate in the polymer. These equations were originally developed to 
compute the high temperature, inelastic deformation of metals. However, variations of these 
equations and the similar set of equations developed by Bodner [61] have been used to model the 
inelastic response of polymers [21,62], A tensorial state variable was used to represent an 
“internal stress’ 7 which represented the resistance to inelastic deformation. As discussed earlier, 
in metals this “internal stress’ 7 was defined as the “back stress”. In the analysis of polymers, the 
“internal stress” represented the resistance to molecular flow. The tensorial state variable is 
orientation dependent. The state variable was assumed to be equal to zero when the material is in 
its virgin state, and evolved towards a maximum value at saturation. As a reminder, the 
“saturation stress” in metals and the “yield stress” in polymers have similar definitions. In this 
report, when the terms “saturation” or “saturation stress” are used in relation to polymers, the 
“yield stress” is the physical mechanism that is being referred to. However, in order to maintain 
consistent terminology with the original development of the equations, the term “saturation” will 
be used. 

The inelastic strain rate was defined as being proportional to the exponential of the 
overstress, the difference between the applied stress and the internal stress. This type of 
relationship was similar in form to the molecular based constitutive equations for polymers 
discussed by Ward [13]. In these equations, the inelastic strain rate (equal to the total strain rate 
at yield) was defined as being proportional to the exponential of the yield stress. 

In the uniaxial simplification of the Ramaswamy-Stouffer constitutive model, the inelastic 
strain rate was defined by the following equation: 


£' = A> ex P| 
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where e' was the inelastic strain rate, a was the stress, and Q was the state variable (“internal 
stress”) which represented the resistance to molecular flow. Material constants included D 0 , a 
scale factor which represented the maximum inelastic strain rate; n, a variable which controlled 
the rate dependence of the deformation response; and Zo, which represented the isotropic, initial 
“hardness” of the material before any load was applied. The value of the state variable Q was 
assumed to be zero when the material was in its virgin state. While not apparent from the 
uniaxial form of the flow equation, Q was a tensorial, not a scalar, variable. 

The evolution of the state variable Q was described by the following expression, 

n = qn m £'- q n |e'| 0.3) 


where Q was the state variable rate,e' was the inelastic strain rate, and Q was the current value 
of the state variable. Material constants included q, which was the “hardening” rate, and O m , 
which was the maximum value of the “internal stress” at saturation. This equation is slightly 
different from the evolution law developed by Ramaswamy and Stouffer [24], The original 
equation included an additional stress rate term that was not used here. The stress rate term was 
included in order to provide for additional hardening at low stress levels. This additional 
hardening was found not to be required for the materials analyzed in this study. For tensile 
loading, where the absolute value of the inelastic strain rate equaled the inelastic strain rate, 
Equation (3.3) was integrated to obtain the following form: 

& = &„, -O m exp(-4£ / ) (3.4) 

where e 1 was the inelastic strain and all other parameters were as defined for Equation (3.3). 

3.3 Material Constant Determination 


A summary of the procedure for determining the material constants for the Ramaswamy- 
Stouffer model will be described here. Detailed discussions of the methods for finding the 
constants can be found in [24,61,63], D 0 was assumed to be equal to a value 10 4 times the 
maximum applied total strain rate, and was considered to be the limiting value of the inelastic 
strain rate. Future investigations may be conducted to investigate whether a relationship between 
D 0 and the shear wave speed can be determined. 

To determine the values of n, Zo and f2 m , the following procedure was utilized. First, the 
natural logarithm of both sides of Equation (3.2) was taken. The values of the inelastic strain 
rate, stress, and state variable Q at “saturation” were substituted into the resulting expression. 
The following equation was obtained; 


In 





V3e 0 

2D 0 
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J 


= In ln(Z 0 ) - 2« ln(<j t - Q m ) 


(3.5) 


where g s equaled the “saturation” stress, e, was the constant applied total strain rate, and the 
remaining terms were as defined in Equations (3.2) and (3.3). 
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The required constants were determined from a set of tensile curves obtained from constant 
strain rate tests. Each curve in this set was obtained at a different constant strain rate. Data pairs 
of the total strain rate and saturation stress values from each curve were taken. Values for Q m 
were estimated for the material, with initial estimates ranging from 50% to 75% of the highest 
saturation stress found to work well. These estimates were similar to the values used for large- 
grain metals. For each strain rate, the data values were substituted into Equation (3.5), and 
represented a point on a master curve. The number of points in the master curve equaled the 
number of strain rates at which tensile tests were conducted. A least squares regression analysis 
was then performed on the master curve. As suggested by Equation (3.5), the slope of the best fit 
line was equal to -2*n. The intercept of the best fit line was equal to 2*n*ln(Z 0 ). The value for 
Q,„ was then adjusted until an optimal fit to the data was obtained. 

To determine the value for q, Equation (3.4) was utilized. At saturation, the value of the 
internal stress was assumed to approach the maximum value, resulting in the exponential term 
approaching zero. Assuming that saturation occurred when the following condition was satisfied: 

exp(-^e/ ) = 0.01 (3.6) 

the equation was solved for q, where 8 S ' was the inelastic strain at saturation. The inelastic strain 
at saturation was estimated by determining the total strain at saturation from a constant strain rate 
tensile curve, and subtracting the elastic strain. The elastic strain was computed by dividing the 
saturation stress by the elastic modulus, as suggested by Equation (3.1). The computed inelastic 
strain was substituted into Equation (3.6), which was solved for q. If the inelastic strain at 
saturation was found to vary with strain rate, the parameter q was computed at each strain rate 
and regression techniques utilized to determine an expression for the variation of q. 

An important point to note is that the meaning and ranges of the material constants when 
these equations were used to model polymer deformation were not the same as when they were 
applied to metals. For example, when analyzing metals using the Ramaswamy-Stouffer model, a 
value of n greater than or equal to 3.0 indicated that the deformation response was relatively rate 
insensitive [24], However, when the equations were applied to polymers, this rule of thumb did 
not necessarily apply due to the differences in deformation mechanisms and strain rate sensitivity 
between the two types of materials. Similarly, other rules of thumb utilized in estimating and 
interpreting the material constants for metals most likely would not apply to polymers. 

3.4 Three-Dimensional Extension of Constitutive Equations 

This section describes the extension of the one-dimensional constitutive equations (Equations 
(3.2)-(3.4)) to three dimensions. To account for the effect of hydrostatic stresses on the response, 
the effective stress term in the flow equation was appropriately modified. A discussion of the 
nature of the tensorial state variable is presented, and the procedure for determining the material 
constants for the three-dimensional version of the equations is discussed. 
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3.4.1 Original Flow Equation 

The three-dimensional extension of the Ramaswamy-Stouffer flow equation given in 
Equation (3.2) was as follows [24]: 


K - D „ exp 


( Z 2 v ' 


3K, 


2 ) 



(3.7) 


where Sjj was the deviatoric stress component, Qy was the component of the state variable, e! 

was the component of inelastic strain, n and Z Q were as defined in Equation (3.2), and K 2 was 
defined as follows: 



where repeated indices indicated summation using the standard indicial notation definitions [24], 
3.4.2 Modified Flow Equation with Shear Correction Factor 

When the flow law described in Equations (3.7) and (3.8) was implemented into the 
composite micromechanics method described later in this report, the nonlinear deformation 
response and component stresses for laminates with shear dominated fiber orientation angles 
were predicted incorrectly. These results indicated that the nonlinear shear stresses in the matrix 
constituent were not being computed properly. For this work, the cause of the discrepancy was 
assumed to be due to the effects of the hydrostatic stresses on the inelastic response. As 
described earlier in this report, hydrostatic stresses have been found to affect the yield behavior 
of polymers. Furthermore, earlier work has indicated that the effective stress terms may need to 
be modified for the multiaxial analysis of polymers [18]. 

The effect of the hydrostatic stresses on the inelastic response of the polymer was accounted 
for by modifying the effective stress term K 2 in the flow law (Equations (3.7) and (3.8)). 
Specifically, since the shear response of the polymer required modification, the shear terms in 
the effective stress were adjusted to account for the effects of hydrostatic stresses. Equation (3.8) 
was rewritten as follows: 


K 2 — — [AT,, + K 22 + K 3? + 2(AT i2 + K n + K 22 )] (3.9) 

The normal terms (1 1,22,33) in this expression maintained their original definition as suggested 
by Equation (3.8) as follows: 


*n -(Sn-GnXSn-fl,,) (3.10) 

^22 =( S 22 - Q 22 X ^22 - Q 2 2 ) (3.11) 
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(3.12) 


However, the shear terms in the effective stress definition were modified as follows: 

k 12 =a(s l2 -n l2 Xs l2 -n ]2 ) 
A' 1 ,=flf(5 13 -n i ,X5 13 -« 13 ) 

K 2} =a{s 2? -Q 2} \S 22 -Q 23 ) 


where: 


a = 
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(3.17) 


J,=jS„S;, (3.18) 

The primary modification to these equations was the multiplication of the shear terms in the 
effective stress by the parameter a. In Equation (3.16), c m was the mean stress, J 2 was the second 
invariant of the deviatoric stress tensor, and P was a rate independent material constant. While 
this formulation was somewhat phenomenological in nature, it was based on actual observed 
physical mechanisms. When the parameter p is set equal to zero (0), the value of a in Equation 
(3.16) is equal to one (1), and Equation (3.9) is equivalent to Equation (3.8). Therefore, the 
modification to the constitutive equations was implemented through the use of the correlation 
coefficient a. 

Since only uniaxial tensile data were available for the polymers considered in this study, the 
value of the parameter p was determined empirically by fitting composite data with shear 
dominated fiber orientation angles, such as [10°] or [15°]. Analyses of [45°] laminates, in which 
the magnitudes of the normal and shear stresses were relatively close, were then conducted in 
order to verify that the determined constant value was reasonable. Ideally, the polymer model 
would be characterized by using a combination of tension, torsion, and tension-torsion tests done 
on the bulk polymer. Since composite data were utilized to characterize the model, simplified, 
consistent techniques to characterize the polymer using bulk polymer data have not yet been 
determined. Future work will involve obtaining multiaxial test data, and determining structured, 
consistent means of determining the effects of hydrostatic stresses on the polymer deformation. 
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3.4.3 Three-Dimensional Extension of Internal Stress Evolution Law 


The internal stress rate, Q tl , was defined by the following relation: 


a/ 






where e' was the effective inelastic strain rate, defined as follows: 


(3.19) 




(3.20) 


and repeated indices again indicated summation using the standard indicial notation definitions 
[24], The remainder of the terms were as defined in Equations (3.2) and (3.3). 

3.4.4 Tensorial Definition of Internal Stress State Variable 


A key difference between the three dimensional Ramaswamy-Stouffer equations, and other 
equations that use tensorial state variables, such as the Walker model [64], lies in the definition 
of the tensorial state variable. In the Ramaswamy-Stouffer model, the back stress (or internal 
stress for polymers) was defined in the same manner as the stress deviator. Under uniaxial 
loading, the tensorial state variable tensor had the following format: 





0 


0 




0 



(3.21) 


where Q represented the uniaxial value of the state variable. In the Walker model, the state 
variable tensor under uniaxial loading was defined as follows. This definition was similar to the 
definition of the plastic strain tensor in standard plasticity theory [24]: 


[ 0*1 = 


Q 

0 


0 


0 



0 


0 
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(3.22) 


While both tensors were deviatoric tensors, the varying definition of the state variable 
resulted in two definitions of the K 2 term. Using the definition of Qj, specified in Equation 
(3.22), K 2 would be defined as follows: 
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(3.23) 


which was what was used in the Walker equation. Note that these expressions were based on the 
original definitions of the effective stress K?, without the shear correction factor discussed above. 

3.4.5 Material Constants for Three-Dimensional Extension of Constitutive Equations 

As discussed by Stouffer and Dame [24], the values of the material constants for the three- 
dimensional formulation of the equations were identical to those obtained using the uniaxial 
representation of the equations. However, since the shear correction factor was not utilized for 
uniaxial tensile loads, the value of the material constant [3 associated with this term would have 
to be determined through multiaxial tests of the polymer. Alternatively, as discussed above (3 
could be determined empirically through examination of tensile curves obtained from off-axis 
loading of composite laminates containing the polymer under consideration. 

3.5 Numerical Implementation of Constitutive Equations 

To test and correlate the constitutive equations, a stand-alone computer code was developed. 
To integrate the flow and evolution laws, the standard fourth order Runge-Kutta explicit 
integration routine was used [65], For this class of equations, implicit integration routines have 
often been used because of their inherent numerical stability [24]. However, in this case, the 
equations were ultimately implemented into a transient dynamic finite element code, which used 
explicit integration schemes. Therefore, in developing the stand-alone computer code, an explicit 
integration scheme was used in order to facilitate the eventual finite element implementation. 
The Runge-Kutta method was employed for this study due to its simplicity and ease of 
implementation. Future efforts might include investigating more robust numerical techniques 
such as semi-implicit algorithms, which provide the stability of implicit methods while still 
maintaining the appearance of an explicit technique. 

To compute the value of a set of variables y n at time step t+At, where t is the current time and 
At is the time increment, the following equation was used: 


V„ (t + At) = y (t ) + — {k ] + 2 k-, + 2 k 4 +k 4 ) 
6 

*i = At * y' (t, y „ ) 


k , = A t* y' 


' 1 . 1 , A 
r + -A/,y„ +-k { 


k, = At* y' 


k 4 = At* y' 


t +— At, y„ + — k 1 


1 


1. ^ 


r + -Ar, v„ +-k. 


where y' was the time derivative of variable y n . 


(3.24) 

3.25) 

(3.26) 

(3.27) 

(3.28) 
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For the stand-alone code developed to correlate and test the constitutive equations, strain 
controlled loading was assumed. This condition was applied for two reasons. First, in a finite 
element application, strains and/or strain increments were passed from the main code into a user 
developed material subroutine. The routine then computed the stresses based on the supplied 
strains. Since the equations utilized here were implemented within a transient dynamic finite 
element code, similar conditions were assumed here to assure compatibility. Second, the tensile 
tests that were used to correlate the model were conducted at constant strain rate. Utilizing strain 
controlled loading in the stand-alone computer code simplified the simulation of these tests. 

To determine the value of the total strain, inelastic strain, and internal stress at time t+At, the 
following algorithm was utilized for each step of the Runge-Kutta integration. The strains at time 
t (or strain estimate at time t+0.5At) were passed into the integration routine. The stresses were 
computed using the elastic constants and the current value of the inelastic strains. The effective 
stress Ko was then computed using Equations (3.9)-(3.18). The inelastic strain rate was computed 
using Equation (3.7), and the internal stress rate was computed using Equations (3.19) and 
(3.20). The elastic Poisson’s ratio and the inelastic strain rates were then used to compute the 
total strain rates. The total strain, inelastic strain and internal stress were computed by integrating 
the corresponding rate equations using Equations (3.24)-(3.28). 

3.6 Model Correlation Analyses 


In this section, the correlation and characterization of the constitutive equations for two 
materials, Fiberite 977-2 and PEEK, is described. 

3.6.1 Fiberite 977-2 Toughened Epoxy 

Consider first Fiberite 977-2, a toughened epoxy. Since toughening the polymer makes it 
more ductile and damage resistant [66] than traditional epoxy resins, a material of this type 
would be more likely to be used in applications where impact resistance is required. 

Tensile tests at low strain rates, ranging from lxlO' 4 /sec to 0.1 /sec, were conducted on the 
Fiberite 977-2 material by Cincinnati Testing Labs, Inc. [67]. Although the ultimate application 
of this research is the study of high strain rate impact problems, good high strain rate tensile data 
were not obtainable at this time. However, the polymer under consideration was rate dependent 
even at relatively low strain rates. Therefore, the ability of the constitutive equations to capture 
rate dependent, inelastic deformation was capable of being examined even by using low strain 
rate data. Furthermore, even though these initial correlation studies were conducted at low strain 
rates, there is no reason, once appropriate data are obtained, that the methodology cannot be 
extended to high strain rate applications. 

Axial tensile tests at constant strain rate were conducted to obtain the test data. Engineering 
stress and engineering strain were measured due to the small strain assumptions that were made. 
Stress-strain curves obtained at constant strain rates of lxlO 4 /sec, 0.01 /sec and 0.1 /sec are 
shown in Figure 3.1. As can be seen in the figure, the tensile response was rate dependent. An 
important point to note is that the unusual results observed at low strains for the results obtained 
at strain rates of 0.01 /sec and 0.1 /sec were most likely due to difficulties that were encountered 
in ramping up to the desired strain rate during the experimental tests. 
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To obtain the tensile material constants, the stress-strain curves had to be extrapolated since 
the tensile specimens failed before “saturation” occurred. To obtain an estimate of the saturation 
stress and strain, the curves for the strain rates of lxlO’ 4 /sec and 0.1 /sec were extrapolated using 
a quadratic curve fit. Since the specimens failed in tension before the saturation level was 
reached, it was assumed that matching the actual saturation stress would not be critical. By 
extrapolating, estimates of the saturation stress and strain were obtained and are shown in Table 
3.1. The elastic and inelastic material constants for this material are shown in Table 3.2. Note 
that the constant (3, related to the shear correction factor, was not obtained using the uniaxial 
tensile data, but is included here for completeness. The procedures used to determine the value of 
this constant will be discussed in a later section of this report. 

Tensile stress-strain curves were computed for constant strain rates of lxlO’ 4 /sec, 0.01 /sec, 
and 0.1 /sec using Equations (3.7)-(3.20). The computed results, along with the experimental 
data, are shown in Figures 3. 2-3.4. As can be seen in the figures, the computed results correlated 
reasonably well with the experimental values. In particular, the nonlinearity of the stress-strain 
curves was captured for all three strain rates. Any discrepancies between the experimental and 
computed results were most likely due to the approximations that were made in estimating the 
saturation stress and strain. The constants used in the constitutive equations were strongly 
dependent on these parameters. Therefore, any inaccuracies in these values could significantly 
affect the computed results. Furthermore, for the experiments conducted at strain rates of 0.01 
/sec and 0.1 /sec, as mentioned earlier difficulties were encountered during the tests in ramping 
up to the desired strain rate. These difficulties could be the cause of the discrepancies between 
the experimental and computed results observed at the lower strain levels. 

3.6.2 PEEK Thermoplastic 

To further examine the capabilities of the constitutive equations, a PEEK 
(polyetheretherketone) thermoplastic was characterized and modeled. Tensile stress strain curves 
were obtained by Bordonaro and Krempl [68] at constant strain rates ranging from 1x10 6 /sec to 
lxlO’ 3 /sec. The experimental results are shown in Figure 3.5. As can be seen in the figure, the 
tensile curves for this material exhibited a distinct “saturation”, or flattening out, unlike the 
tensile curves for the Fiberite 977-2 epoxy. This result was expected since a thermoplastic like 
PEEK is more ductile than the toughened epoxy Fiberite 977-2. Due to the fact that the tensile 
curves flatten out, the inelastic material constants were determined directly, with no 
approximation (beyond those standard to the constitutive equations) required. The inelastic 
material constants for this polymer were determined using the procedures described earlier, and 
are shown in Table 3.2 along with the elastic modulus and Poisson’s ratio. 

Tensile stress-strain curves were computed under strain controlled loading at constant strain 
rates of lxlO’ 6 /sec, lxlO’ 4 /sec and lxlO' 3 /sec. The computed curves and the corresponding 
experimental results are shown in Figures 3. 6-3. 8. As can be seen from the figures, the computed 
values correlated with the experimental results extremely well. Furthermore, the constitutive 
equations appeared to capture the tensile response of this polymer more accurately than was seen 
for the Fiberite 977-2. The quality of the simulations for the PEEK was most likely due to the 
fact that the measured stress-strain curves for this material displayed a distinct saturation, 
allowing for a more accurate determination of the material constants. These results indicated that 
if the polymer under consideration was appropriately characterized, the constitutive equations 
described here could do a good job in computing the polymer deformation response. 


NASA/TM — 1 999-209768 


21 



3.7 Summary 


Rate dependent, inelastic constitutive equations based on the state variable method have been 
used to model the deformation response of ductile polymers under small strain conditions. Two 
representative polymers were characterized, and their uniaxial tensile deformation response was 
computed using the constitutive equations. The computed results correlated reasonably well to 
the experimental values for both polymers, indicating that the equations adequately captured the 
material response. The constitutive model can thus be used in micromechanics equations to 
predict the rate dependent, nonlinear deformation response of polymer matrix composites. 


Table 3.1 

Extrapolated Saturation Stress and Saturation Strain Values for Fiberite 977-2 


Strain Rate (/sec) 

Saturation Stress (MPa) 

Saturation Strain 

lxlO' 4 

97 

0.055 

0.1 

110 

0.053 


Table 3.2 

Materia] Properties for Fiberite 977-2 and PEEK 



E 

(GPa) 

V 

D 0 

(1/sec) 

N 

Zo 

(MPa) 

q 

(MPa) 

P 

977-2 

3.65 

0.40 

1E+04 

0.50 

1030 

100 

69 

1.2 

PEEK 

4.00 

0.40 

1E+04 

0.70 

630 

310 

52 

0.45 
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Figure 3.1: Tensile Curves for Fiberite 977-2 Toughened Epoxy 
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Figure 3.2: Model Correlation for 977-2 Resin at Strain Rate of lxlO' 4 /sec 
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Figure 3.5: Tensile Curves for PEEK 
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CHAPTER 4 

COMPOSITE MICROMECHANICAL MODEL 


The implementation of the polymer constitutive equations described above into a mechanics 
of materials based micromechanics model will now be described. In this manner, the nonlinearity 
and rate dependence of the deformation response of a polymer matrix composite were accounted 
for by computing the inelastic response of the polymer matrix. In this chapter, the assumptions of 
the micromechanics model will be presented, along with an overview ol the method. The 
equations used to compute local and effective stresses as a function of total strains will be 
derived, and the equations used to compute the effective inelastic strains will be given. The 
numerical techniques used to implement the micromechanics equations into a stand-alone 
computer code will be discussed. Finally, verification studies conducted using two representative 
polymer matrix composites will be described. 

4.1 Model Assumptions 

Micromechanics techniques were used to predict the effective properties and deformation 
response of the polymer matrix composites examined in this study. As mentioned previously, in 
micromechanics methods the effective properties of a composite material unit cell are predicted 
based on the properties of the individual constituents. As was discussed earlier in this report, 
micromechanics techniques were utilized for this study since the inelastic, rate dependent 
deformation response of the polymer was explicitly modeled using the constitutive equations 
described in the previous chapter. Furthermore, a fairly simple mechanics of materials based 
composite micromechanics model was utilized in order to facilitate implementation of the model 
within a transient dynamic finite element computer code. 

For this study, the composite unit cell was defined as consisting of a single continuous fiber 
and its surrounding matrix. Only laminated composites were analyzed; woven composites were 
not considered at the present time. The matrix constituent and the composite as a whole were 
assumed to have a sufficient degree of ductility for nontrivial inelastic strains to be present. 

The composites were assumed to have a periodic, square, fiber packing arrangement, with 
perfect bonding between the fiber and the matrix. These assumptions are common in the 
micromechanical analysis of composite materials [33-36]. While actual composites often have 
more complicated fiber architectures [69], for this study the fiber packing arrangement utilized 
was chosen in order to simplify the development of the micromechanics equations and to 
minimize the computational effort required. If future analytical results indicate that fiber packing 
plays a significant role in the application under consideration, modifications to the 
micromechanical models could be made. Alternatively, selected detailed finite element analyses 
might be performed in order to quantify the fiber packing effects. The assumption of perfect 
bonding, a common assumption for polymer matrix composites, was also made in order to 
simplify the development of the micromechanics equations. If fiber/matrix debonding turns out 
to play a significant role in the strength and failure analyses of the materials under consideration, 
once again appropriate modifications could be made to the equations. 

Only unidirectional composites at various fiber orientation angles were analyzed with the 
micromechanics equations presented in this study. The finite element method could be used to 
model laminated composites with varying fiber orientation through the thickness. In such an 
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analysis, two approaches would be possible. First, a layer of elements could be used to model a 
single ply of the composite at a specified orientation angle. Multiple layers of elements would 
then be used to simulate the composite laminate. Alternatively, if the finite element code under 
consideration permitted laminated shell elements, the micromechanics equations presented in 
this section would be used to calculate the stresses for a single layer of the laminated element. 

As discussed earlier, some efforts have been made by previous researchers to apply equations 
of state on the micromechanical level in modeling the high strain rate response of polymer 
matrix composites [44]. Equations of state were used to model the effects of changing density on 
the hydrostatic stresses in the material. These equations are usually only required for very high 
strain rate loading conditions. The strain rates encountered in fan containment problems should 
be low enough that equation of state considerations on the micromechanical level will not be 
required. 

The deformation response of the polymer matrix was simulated using the modified 
Ramaswamy-Stouffer constitutive equations described above. The fibers of the composite were 
assumed to be linear elastic, with rate independent properties. Temperature effects were not 
considered, and small strain conditions were assumed. These assumptions were also applied in 
the development of the matrix constitutive equations. 

4.2 Overview of Micromechanics Method 


The micromechanics method utilized in this study was based on a method proposed by Sun 
and Chen [70]. In this approach, the composite unit cell was broken up into three subcells. One 
subcell represented the fiber while the remaining two subcells represented the matrix. This 
approach was similar to the Method of Cells approach utilized by Aboudi [32]. However, in the 
Method of Cells a displacement field was assumed for each subcell. The effective displacements 
and stresses were then determined using the equations of continuity and equilibrium. As part of 
this process, the stresses in the individual subcells were also computed. In the Sun and Chen 
approach, on the other hand, uniform stress and uniform strain assumptions, in combination with 
the material constitutive equations, were utilized to solve for the stresses and strains for each 
subcell and for the overall composite. Furthermore, in the Sun and Chen model out of plane 
stresses were neglected (plane stress), and classical plasticity theory was used to account for any 
inelastic strains which might be present. In the Sun and Chen model, stresses were assumed to be 
given, and the strains were determined by the constitutive model. 

Robertson and Mall revised and expanded the Sun and Chen model [71-73]. In this approach, 
the plane stress assumption was removed, and the full three-dimensional stresses and strains 
were computed for each subcell and for the overall composite. Since the model was fully three- 
dimensional, four subcells were used to represent the unit cell. One subcell represented the fiber, 
and the remaining three subcells represented the surrounding matrix material. The Robertson and 
Mall model applied unified state variable constitutive equations to compute the inelastic strains 
in the matrix material. The micromechanics equations were formulated assuming stress 
controlled loading, in which the subcell stresses and strains were computed based on a defined 
effective stress condition. Since Robertson and Mall concentrated on analyzing metal matrix 
composites, where fiber/matrix debonding is significant, the equations were also modified to 
allow for the presence of a weak fiber/matrix interface. 
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Pindera and Bednarcyk [74] utilized a similar approach in reformulating the Generalized 
Method of Cells [75], The Generalized Method of Cells was a reformulation of the original 
Method of Cells [32] that allowed for an unrestricted number of subcells in the unit cell. In this 
latest reformulation by Pindera and Bednarcyk, the subcell stresses were computed based on the 
macroscopic strains. The reformulation significantly improved the computational efficiency of 
the method. 

The micromechanics model developed for this study was similar to the method utilized by 
Robertson and Mall. Uniform stress and uniform strain assumptions were applied to a four 
subcell unit cell, which is displayed in Figure 4.1. In the unit cell, subcell “Af” represented the 
fiber and subcells “Am”, “Bl”, and “B2” were composed of matrix material. The material axis 
system was as shown in the figure. The “1” coordinate direction was along the fiber direction, 
the “2” coordinate direction was perpendicular to the fiber in the plane of the composite, and the 
“3” coordinate direction was perpendicular to the fiber in the out of plane direction. The fiber 
was idealized as having a square shape, with the side length equal to the square root of the fiber 
volume fraction (k t ). Assuming a square fiber shape would result in an incorrect prediction of the 
interfacial stresses. However, due to the perfect bonding assumption, as well as the expected 
failure modes in the chosen application, the accurate prediction of interfacial stresses was 
assumed not to be critical. 

The full three-dimensional stresses and strains were computed for each subcell and for the 
unit cell as a whole. By removing the plane stress assumption, thick composites could be 
analyzed. Furthermore, through the thickness stresses could be more accurately computed, which 
would most likely be important in modeling high strain rate impact normal to the plane of the 
laminate. 

In the equations developed here, the subcell stresses and equivalent unit cell stresses were 
computed based on the applied macroscopic strains. The loading condition was the primary 
difference between this method and the Robertson and Mall technique. Utilization of strain 
controlled loading simplified the implementation of this model into a finite element code. In a 
user defined material subroutine in a finite element code, strains were passed into the routine, 
and stresses were computed and passed back to the calling routines. Furthermore, the equations 
presented here were to some extent a specialization of the reformulation of the Generalized 
Method of Cells (GMC) described above. Specifically, if the GMC equations were explicitly 
solved for a four subcell model, equations very similar to those presented in this report would 
result. However, while GMC was designed to allow for a variable number of subcells, the 
equations given in this report were limited to a four subcell unit cell model, which simplified the 
equations significantly. 

4.3 Derivation of Micromechanics Equations 

The unit cell used in the development of the micromechanics equations is shown in Figure 
4.1. The bottom layer of subcells, with subcells “Af” and “Am”, was referred to as Row 1 ( R 1 ) . 
The top layer of subcells, with subcells “Bl” and “B2”, was referred to as Row 2 (R2). Column 1 
(Cl) was defined as consisting of subcells “Af’ and “Bl”, and Column 2 (C2) was defined as 
consisting of subcells “Am” and “B2”. The subscript “f’ was used to denote fiber related 
properties, and the subscript “m” was used to denote matrix related properties. Subscripts “Af’, 
“Am”, “Bl” and “B2” were used to denote stresses and strains of the individual subcells. 
Subscripts “Rl”, “R2”, “Cl”, and “C2” were used to denote stresses and strains in the 


NASA/TM — 1 999-209768 


29 



corresponding regions as defined above. Stresses and strains with no region identifying subscript 
were assumed to represent the total effective stresses and strains for the unit cell. A superscript 
“I” was used to denote inelastic strains. The subscripts “11”, “22” and “33” were used to define 
normal stresses, strains and material properties, with the coordinate directions as defined in 
Figure 1. The subscripts “12”, “13” and “23” were used to define shear stresses, strains, and 
material properties. 

The symbol “E” represented the elastic modulus, the symbol “G” represented the shear 
modulus, and the symbol “v” represented the Poisson’s ratio. Subscripts were attached to these 
terms as noted above. The symbol “Gi” represented stress tensor components, the symbol “£y” 
represented strain tensor components, and the symbol “y,” represented engineering shear strain 
components, all assigned in a Cartesian frame of reference. The symbol “kf” represented the fiber 
volume ratio of the composite. 

The stress and strain in each subcell were assumed to be the effective stress and strain, equal 
to the average stress or strain over the volume of the subcell, and were assumed to be uniform 
over the volume of the subcell. The effective stress and strain in Row 1, Row 2, Column 1 and 
Column 2 were defined as the volume average of the stresses and strains in the component 
subcells. The effective stress and strain in the unit cell were defined as the volume average of the 
stresses and strains in Row 1 and Row 2 (or Column 1 and Column 2). To determine the volume 
average, a weighted sum was computed where the value (stress or strain) in each subcell or 
combination of subcells was weighted by the volume ratio of the subcell or combination of 
subcells. 

The fibers were assumed to be transversely isotropic. The components in the transversely 
isotropic compliance matrix (the inverse of the stiffness matrix) for the fiber were defined as 
follows. Note that the symbol Sjjf, which was used to denote the terms in the compliance matrix 
for the fibers, should not be confused with the symbol Sjj, which was used in the previous 
chapter to represent the components of deviatoric stress in the polymer constitutive equations. 


S = — 
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' 22 / 
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where Enf represented the longitudinal elastic modulus of the fiber (along the 1 direction axis in 
Figure 4.1), Ei 2 f represented the transverse elastic modulus of the fiber, Vi 2 f represented the axial 
Poisson’s ratio of the fiber, V 23 f represented the transverse Poisson’s ratio of the fiber, G^f 
represented the in-plane shear modulus of the fiber, and G 2 3 f represented the transverse shear 
modulus of the fiber. 

The matrix was assumed to be an isotropic material, with compliance matrix terms Sij m 
defined as follows. 


5 


llm 






66 m 



(4.2) 


where E m represented the elastic modulus of the matrix, v m represented the Poisson’s ratio of the 
matrix, and G m represented the shear modulus of the matrix. 
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The transversely isotropic compliance matrix was used to relate the strains to the stresses, 
using the following relations. Again, note that in these equations Sy represented the components 
of the compliance matrix, not the components of the deviatoric stress tensor as in the previous 
chapter. 
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The addition of the inelastic strain components to the standard transversely isotropic elastic 
constitutive law was how inelasticity was incorporated into the constitutive relations. For the 
fiber, which was assumed to be linear elastic, these components were neglected. For the matrix 
material, which was assumed to be isotropic, Si? was set equal to S|i, S 22 was set equal to Sn, 
and S 44 was set equal to S& 6 . 

The effective total strains in the composite unit cell were assumed to be known in the 
micromechanics equations. Furthermore, the inelastic strains in each subcell were assumed to be 
given. Also, the MATHCAD software package [76] was utilized to assist in carrying out the 
algebraic computations presented in this derivation. 

4.3.1 Normal Stresses and Strains 


For the normal stresses and strains ( 1 1, 22 and 33), the following uniform stress and uniform 
strain assumptions were made: 


In the fiber direction: 
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Normal to the fiber, in the plane of the ply: 
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Normal to the fiber, normal to the plane of the ply: 


^33,4/ ^3301 — ^*330 

(4.10) 

^*3 3,4/// = ^33fi2 = ^33C2 


^33C1 = ^ 33C2 ~ ^33 

(4.11) 

The effective stresses and strains in Row 1, Row 2, Column 1 and Column 2, as well as for 

the composite unit cell, were computed using volume averaging, 
expressions: 

yielding the following 

£ 22R\ - * ^22 At + (1 “ 2 Am 

(4.12) 

^22R2 ~ #7 ^22B] 0 — \j ^ f ) * ^22B2 

(4.13) 

^33C1 = -\]k f ' ^33 Af — ^ f ) * ^3351 

(4.14) 

^33C2 “ ^ ^33Am + “ yj ^ f ) * ^33B2 

(4.15) 

~ ' &\ \Af V^7 ^ * ^UAm 

(4.16) 

^ll/?2 — * ^llfll + (1 — y[^"f } ' ®\)B2 

(4.17) 

~ yjk t ‘ r\ + (1 ” ijkf ) * &}]R2 

(4.18) 

® 22 7^/ ' ^ 22R 1 (1 “* ijk f ) * tT 22R2 

(4.19) 

^33 — 7^7 ' ^?3C1 (1 ~ #7 ^ ' ^33C2 

(4.20) 

The constitutive relations for the fiber and matrix were defined as follows, using Equation (4.3): 
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By solving Equations (4.21) and (4.24) for each subcell, and by utilizing the appropriate 
uniform stress and uniform strain assumptions, the following expressions were obtained. 
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Equations (4.27)-(4.30) were then substituted into Equations (4.22), (4.23), (4.25), and (4.26) 
for each subcell, applying the appropriate uniform stress and uniform strain assumptions. By 
inserting the resulting expressions into Equations (4.12)-(4.15), the following system of 
equations resulted: 
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Equations (4.31)-(4.34), together with Equations (4.27)-(4.30), were solved for the required 
subcell stresses. Equations (4.8), (4.9), and (4.16)-(4.20) were then used to compute the effective 
stress state in the unit cell. 

4.3.2 In-Plane Shear Stresses and Strains 

For the in-plane shear (1-2 direction) stresses and strains, the following uniform stress and 
uniform strain assumptions were made: 


In the plane of the ply: 


Y\2R\ — Y \2R2 ~ Y 12 

(4.35) 

®\2Af ~ ®l2Aw = &12R] 
&\2B\ = ®\2B2 = ®]2R2 

(4.36) 

By applying volume averaging, the effective in-plane shear stresses and strains for Row 1, Row 
2, and the composite unit cell were defined as follows: 

YllRl - * Y\2Af + (1 #7 ^ * Y\2Am 

Y\2R2 ~ yjk f * 7l2Bl + (1 “ -yjk f ) * 7j 2 fi2 

(4.37) 

^12 2R\ + (1 f ) * ^*12/? 2 

(4.38) 
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The constitutive relations for the fiber and matrix were defined by the following expressions, 
using Equation (4.4): 


12/ “ ^66; ' ^12/ 

(4.39) 

— S * cr 

°66 /H c l2/» 

(4.40) 


By substituting Equations (4.39) and (4.40) into Equation (4.37), and using the appropriate 
uniform stress and uniform strain assumptions, the following expressions were obtained, from 
which the subcell in-plane shear stresses were computed. 


7 .2 



i 5 66 J*<T, 2 „+2*(l- > /*7)*e.^ 


(4.41) 


7.2 


— S(,bm * ®\2R2 + ~ * bj ^1 " £ .2fil + 2 ' (1 ^Jk f ) £ |2 B2 


(4.42) 


4.3.3 Transverse Shear Stresses and Strains: 1-3 Direction 

The computation of the subcell shear stresses in the 1-3 direction was very similar to the 
computation of the subcell shear stresses in the 1-2 direction. The only difference was that 
Column 1 and Column 2 were used instead of Row 1 and Row 2. The uniform stress and uniform 
strain assumptions shown in Equations (4.35) and (4.36) were transformed as follows: 

In the fiber/normal direction: 

Toci = 7i?c: = 7i3 < 4 - 43 ) 


<^13 At ~ ^130 

= ^13B2 = *^13C2 


(4.44) 


The volume averaged stresses and strains in Column 1 and Column 2 were computed using the 
expressions: 


Toe. -^'•7,^(1 (4 45) 

7 l 3 C 2 = -\jk f y + ( 1 — yj k t ) y lW2 

(T„ = * <T,, C1 <t, 3C2 (4-46) 

By substituting Equations (4.39) and (4.40) (replacing the subscript “12" with “13" as 
indicated by Equation (4.5)) into Equation (4.45), the transverse 1-3 direction shear stresses in 
the individual subcells were computed from the following equations: 

7,., = [-v/^7 * s 66 , + a - fijr * s 66m ] * <7 1? a + 2 * a - V*7> * 4 b. ( 4 - 47 ) 
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(4.48) 


7.3 = S 66,n *Vw+ 2 *4W* + 2 * ( 1 - Jkf ) * *2 

4.3.4 Transverse Shear Stresses and Strains: 2-3 Direction 

To compute the transverse 2-3 direction shear stresses in the subcells of the composite unit 
cell, the following uniform stress assumptions were made: 

In the transverse ply-normal direction: 

^2381 & 21R2 = ® 23 

^23/V = ^23 Am ~ ^2381 
^2381 = ^2382 = ^2382 

The volume averaged shear strains in Row 1, Row 2 and the composite unit cell were then 
defined as follows: 


(4.49) 

(4.50) 


23/?l “ ^ 7:3,4/ + (1 * YllAm 

(4.51) 

23 R2 = f * 72351 + yjk f ) * 7:352 

723 = +(1-V*n *72352 

(4.52) 


The constitutive relations for the fiber and matrix were defined by using Equation (4.6): 

7 23/ — -*44/ ‘ ^23/ (4.53) 

723,„=^66,„ :f: ^23,„+2*4,„ (4.54) 

By substituting Equations (4.53) and (4.54) into Equations (4.51) and (4.52), and by using 
the uniform stress assumptions, the following expression was obtained which was used to 
compute the subcell shear stresses in the transverse 2-3 direction: 

7 2 3 = t^/ * ^44/ + ( J - k f ) * 5^,,, ] * <7,, 

+ 2 * (1 - 7*7) * [^*7 * (4.4,,, + £238. ) + <1 - 7*7) * £ 2382 J (4 ' 55) 

4.4 Effective Inelastic Strains 


In applying a strain controlled load to a composite unit cell, accurate determination of the 
Poisson strains was required. To compute the Poisson strains based on an applied total strain in a 
particular coordinate direction, the effective inelastic strains in each coordinate direction were 
required. To compute the effective inelastic strains for the unit cell, the same uniform stress and 
uniform strain assumptions that were used in the previous section were applied again. 
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The constitutive relations (Equations (4.3)-(4.6)) for the fiber and matrix were modified for 
these calculations. By making these simplifications, the complexity of the mathematics was 
reduced considerably. However, the accuracy of the results was not significantly affected. The 
modified constitutive equations for the normal and in-plane shear stresses were given by the 
following expressions: 


^11 — ^1 1 1 £ 11 ) 

(4.56) 

O -> ■> • — E ->-j (£-)-> 

(4.57) 

o n - G n (y 12 -2*el 2 ) 

(4.58) 


where all symbols were as defined previously. The equations for the transverse shear stresses 
were very similar to Equation (4.58), replacing the subscript “12” with “13” or “23” as 
appropriate. The equation for the stress in the 3-3 coordinate direction was very similar to 
Equations (4.56) and (4.57), replacing the subscript “11” or “22” with “33”. Note, as before, the 
inelastic strain values were only used when the matrix material was considered. Furthermore, for 
the isotropic matrix material, E m replaced E] i and E 22 , and G m replaced G 12 . 

By utilizing the constitutive equations for the fiber and matrix (Equations (4.56)-(4.57)), 
along with the uniform stress and uniform strain assumptions (Equations (4.7)-(4.20)), the 
following expressions were obtained for the effective inelastic normal strains: 

, _ (I- 

£;:= eZfjFW. 


where: 



(4.61) 


The effective inelastic strain in the 3-3 direction was also computed by using Equation (4.60). In 
this case, the inelastic strains in subcells “Am” and “Bl” were switched in the expression. 

By utilizing the constitutive equations for the fiber and matrix (Equation (4.58)), along with 
the uniform stress and uniform strain assumptions (Equations (4.35)-(4.38)), the following 
expression was obtained for the effective inelastic in-plane shear strain: 
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where: 



(4.63) 


The effective inelastic strain in the 1-3 out-of-plane shear direction was computed by using 
Equation (4.62) and switching the inelastic strains in subcells “Am” and “Bl”. The effective 
inelastic strain in the 2-3 out-of-plane shear direction followed automatically from Equation 
(4.55), and is presented without further comment. 

=d-V*7 +4 bi) + (1-V^7 ( 4 -64) 



For the current study, a stand-alone computer code was developed in order to implement and 
test the micromechanics equations. A standard fourth order Runge-Kutta explicit integration 
scheme was again utilized to integrate the rate dependent constitutive equations. The details of 
the method can be found in Section 3.5. 

As mentioned in the development of the micromechanics equations, strain controlled loading 
was assumed in the formulation. In the computer algorithm, strains were specified in a particular 
coordinate direction. To impose the required Poisson and axial-shear coupling strains, effective 
elastic properties for the composite at a specified fiber orientation angle were utilized. First, the 
elastic constants in the material axis system were computed using equations developed by 
Murthy and Chamis [37]. The elastic constants in the structural axis system were then computed 
using standard techniques and equations described in references such as [33], [36] and [47]. The 
material axis system was the coordinate system shown in Figure 4.1. The structural axis system 
was the axis system along which the loads are applied. The material coordinate system was 
obtained by rotating the structural axis system about the “3” coordinate axis by an amount equal 
to the fiber orientation angle. 

For the computer code execution, first the required geometric data (fiber volume ratio and 
fiber orientation angle), constituent properties and load history data were read in from an input 
file. The required elastic constants in both material and structural coordinate systems were 
computed, along with the required tensor transformation matrices. For each time step, the total 
strain rate in the load direction was determined. The Runge-Kutta integration procedure was then 
carried out to compute the total strains in the structural axis system, as well as the inelastic 
strains and internal stresses in each subcell. The total stresses in structural coordinates were 
calculated using the total strains, appropriate tensor transformations, and the micromechanics 
equations. At this point, the code moved on to the next time step. 

The Runge-Kutta integration algorithm involved the computation of several intermediate 
estimates of the total strains, subcell inelastic strains and subcell internal stresses. To calculate 
the intermediate estimates, first the total strain estimate was converted from the structural axis 
system to the material axis system. The stresses in each of the subcells were then determined 
using the micromechanics equations. Using the computed stresses, the inelastic strain rates and 
internal stress rates in each matrix subcell were computed using the polymer constitutive 
equations. The effective inelastic strain rate tensor for the composite unit cell in the material axis 
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system was computed, and the values were then transformed into the structural axis system. 
Using the ply level Poisson’s ratios and axial-shear coupling coefficients, the total strain rate 
tensor in structural coordinates was calculated. The intermediate values required for the Runge- 
Kutta integration routine were then determined. 

4.6 Model Verification Analyses 

To verify the micromechanics equations, a series of analyses were earned out using two 
material systems. Both material systems exhibited a nonlinear deformation response for off-axis 
fiber orientation angles. 

4.6.1 Material Properties 

The first material system examined, supplied by Fiberite, Inc., consisted of carbon IM-7 
fibers in a 977-2 toughened epoxy matrix. Unidirectional laminates with fiber orientations of 
[0°], [10°], [45°], and [90°] were tested. Tensile tests were conducted by Cincinnati Testing Labs 
of Cincinnati, Ohio at a strain rate of 1x10 4 /sec on each of the composites [67], As mentioned 
before in the discussion of the polymer constitutive equations, good high strain rate tensile data 
have not been obtained at this time. However, the equations will be characterized and validated 
at high strain rates once appropriate data are obtained. 

The IM7/977-2 composite had a fiber volume ratio of 0.60. The material properties of the 
IM-7 fibers, as determined from Reference [77], are shown in Table 4.1. The material properties 
of the 977-2 matrix, discussed in Chapter 3, are listed in Table 3.2. As mentioned in Section 
3.4.2, the constant P for the shear correction factor was determined empirically using composite 
data. For the IM7/977-2 system, the constant was determined using the [10°] off-axis composite 
data. This orientation was chosen as it is shear dominated [78]. The value of P was determined to 
be 1.2. 

The second material that was studied consisted of carbon AS-4 fibers in a PEEK 
thermoplastic matrix. Tensile stress-strain curves were obtained by Weeks and Sun [25] for 
unidirectional composites with fiber orientations of [14°], [30°], [45°] and [90°] at a strain rate of 
1x10° /sec, and composites with fiber orientations of [15°], [30°], [45°] and [90°] at a strain rate 
of 0.1 /sec. Only the low strain rate data were examined since the PEEK material was only 
characterized for relatively low strain rates. However, by conducting verification studies using 
this data, the ability of the micromechanics equations to capture the rate dependent deformation 
response of a polymer matrix composite was determined. 

The fiber volume ratio used for the AS4/PEEK material was 0.62 (a typical value for this 
material based on representative manufacturer information). The elastic properties of the AS-4 
fibers, as listed in Reference [79], are shown in Table 4.1. For the PEEK matrix, the material 
properties are shown in Table 3.2. The value of the constant P was determined empirically by 
fitting data from the [14°] and [15°] laminates, which were shear dominated. By examining the 
data, the constant P was determined to be rate independent, and equal to a value of 0.45. Note 
that since there was currently no good explanation for the physical explanation for the physical 
meaning of the parameter p, the significance of the value of this parameter and its variation for 
the two materials examined in this study could not be quantified at this time. 
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4.6.2 Analysis Results 

Analyses were conducted using the micromechanics equations described above. The 
modified Ramaswamy-Stouffer constitutive equations described in Chapter 3 were used to 
compute the rate dependent, inelastic response of the polymer matrix. The predicted results were 
compared to experimentally obtained values. Stress-strain curves for the IM7/977-2 laminates 
are shown in Figures 4. 2-4. 5. Stress-strain curves for the AS4/PEEK composite are shown in 
Figures 4.6-4.9 for a strain rate of 1x10 ' 5 /sec, and in Figures 4.10-13 for a strain rate of 0.1 /sec. 
The calculations for the [10°] IM7/977-2 laminate shown in Figure 4.3, and the [14°] and [15°] 
AS4/PEEK laminates shown in Figures 4.6 and 4.10 were correlations, while the remaining 
calculations were “true” predictions. 

As can be seen in the figures, for both materials, and for both strain rates for the AS4/PEEK 
system, the analytical results matched the experimental values quite well in general for all fiber 
orientation angles examined. For the shear dominated fiber orientation angles ([10°], [14°], and 
[15°]), the inelastic portion of the predicted stress-strain curves were flatter than the experimental 
results. Furthermore, the strain level at which significant inelasticity begins was somewhat 
overpredicted. However, the predicted stresses at the end of the stress-strain curves matched the 
experimental values reasonably well. For predicting ply strength, which will be discussed in the 
next chapter of this report, correctly predicting the stress levels at the end of the stress-strain 
curve was most critical. 

For the AS4/PEEK system, the elastic response of the [30°] and [45°] laminates was slightly 
underpredicted. As the elastic response of the other laminates was predicted reasonably well, 
these results indicated that the transverse Poisson’s ratio or shear moduli used for the AS-4 fiber 
may not have been correct. As these values are often just approximated instead of experimentally 
measured, inaccuracies in these values were not surprising. 

Overall, the comparison between the experimental and predicted values was quite good. In 
particular, the rate dependence and nonlinearity of the deformation response of the composites 
was captured. To emphasize the ability of the analytical model to capture the rate dependence of 
the deformation response, the results for the [30°] AS4/PEEK laminate shown in Figures 4.7 and 
4.11 are superimposed in Figure 4.14. In this figure, the curves labeled “IE-05” are results 
obtained at a strain rate of 1x10'"'' /sec, and the curves labeled 0.1 are results obtained at a strain 
rate of 0. 1 /sec. The results in this figure show that the response of the AS4/PEEK system was 
rate dependent for this fiber orientation. Furthermore, the predicted curves show that the 
micromechanics equations and the polymer constitutive model captured the rate dependence. The 
rate dependence of the composite deformation response should be even more dramatic once high 
strain rate experiments and analyses are considered. 

4.7 Summary 


A mechanics of materials based composite micromechanics model was developed to predict 
the effective properties and deformation response of rate dependent, inelastic, polymer matrix 
composites based on the constituent properties and applied strains. The modified Ramswamy- 
Stouffer equations were used to model the polymer matrix constituent. The effective inelastic 
strains of the composite unit cell have also been computed. The model has been implemented 
into a stand-alone computer code, and was verified for two representative polymer matrix 
composites. The analytical predictions compared favorably to experimentally obtained values. 
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As discussed in the next chapter, the micromechanics model was used in combination with 
selected failure criteria to predict the failure of a composite ply based on the macroscopic 
stresses. 


Table 4. 1 

Material Properties for IM-7 and AS-4 Fibers 



Longitudinal 
Modulus (GPa) 

Transverse 
Modulus (GPa) 

Poisson’s Ratio 

In-Plane Shear 
Modulus (GPa) 

IM-7 

276 

13.8 

0.25 

20.0 

AS-4 

214 

14.0 

0.20 

14.0 



Col. 1 
(Cl) 

Col. 2 
(C2) 


( 1 -sqrt(k f )) 

B1 

B2 

Row 2 (R2) 

sqrt(k f ) 

Af 

Am 

3 

Row 1 (Rl) 

1 

sqrt(kf) 

(l-sqrtfkj)) 



k t =Fiber Volume Fraction 

Af=Fiber Subregion 

Am,B l,B2=Matrix Subregions 

Figure 4. 1 : Geometry and Layout of Composite Unit Cell Model 
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Stress (MPa) Stress (MPa) 



Figure 4.6: Model Correlations for [14°] AS4/PEEK Laminate 
at Strain Rate of lx 10' 5 /sec 


■ Experimental 
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■ Experimental 
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Figure 4.7: Model Predictions for [30°] AS4/PEEK Laminate 
at Strain Rate of lxlO' 5 /sec 
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Stress (MPa) Stress (MPa) 
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Figure 4.8: Model Predictions for [45°] AS4/PEEK Laminate 
at Strain Rate of lxlO' 5 /sec 



0 0.001 0.002 0.003 0.004 0.005 0.006 0,007 

Strain 

Figure 4.9: Model Predictions for [90°] AS4/PEEK Laminate 
at Strain Rate of lx 10 5 /sec 
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Stress (MPa) stress (h 



Figure 4.10: Model Correlations for [15°] AS4/PEEK Laminate 
at Strain Rate of 0. 1 /sec 



Strain 

Figure 4.11: Model Predictions for [30°] AS4/PEEK Laminate 
at Strain Rate of 0. 1 /sec 
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Stress (MPa) Stress (MPa) 
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Figure 4.12: Model Predictions for [45°] AS4/PEEK Laminate 
at Strain Rate of 0. 1 /sec 


60 7 ' 



Strain 

Figure 4.13: Model Predictions for [90°] AS4/PEEK Laminate 
at Strain Rate of 0. 1 /sec 
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Figure 4. 14: Prediction of Strain Rate Dependence of [30°] AS4/PEEK Laminate. 
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CHAPTER 5 

PLY STRENGTH MODEL 


The implementation of a ply failure model into the micromechanics equations developed 
in the last chapter will be presented next. An overview of the methodology used to predict 
ply ultimate strength will be given. The Hashin failure criterion, which was used in this study, 
will be described in detail. The ultimate strengths for the representative polymer matrix 
composites analyzed in the previous chapter will be predicted for a variety of fiber orientations 
and strain rates. 

5.1 Overview 


In order to develop structural level penetration and failure models that can be applied to high 
strain rate impact applications, ply level failure needs to be accurately predicted. As discussed in 
the background section of this report, previous researchers have developed a variety of ply level 
failure models. For this work, ply level failure due to local failure mechanisms was predicted 
based on macroscopic stresses and strengths. In particular, the Hashin failure model [50] was 
utilized in this work. As discussed earlier, some researchers such as Pecknold and Rahman [59] 
have utilized constituent level stresses to predict ply failure based on local failure mechanisms. 
However, for this study, ply level strength data were more readily available and considered to be 
more reliable than constituent level strength data. Nonetheless, in considering the ultimate future 
application of the ply strength models, the ability to account for and predict local failure 
mechanisms in at least an approximate manner was desired. 

Since only ply level failure based on simple tensile tests was considered at this time, property 
degradation models were not utilized. For structural level modeling, the ability to only degrade 
certain material properties based on the local ply failure mechanisms would be desirable to 
provide improved simulation of stress transfer mechanisms. Furthermore, in implementing the 
model into a finite element code, a gradual degradation of the material properties might improve 
the stability of the finite element analysis. Since the failure model utilized in this study predicted 
failure based on approximations of local failure mechanisms, the eventual incorporation of 
property degradation models would be possible. 

5.2 Hashin Failure Criteria 


For this work, the Hashin failure criteria were utilized [50]. These criteria predicted ply level 
failure based on local mechanisms using macroscopic stresses and strengths. These criteria were 
based on stresses and strengths in the principal axes of the material, so appropriate 
transformations [24] were carried out to convert stresses from the structural axis system to the 
principal material axis system. Since the composites analyzed in this study were only subject to 
in plane loading and the out-of-plane stresses were found to be relatively small, the plane stress 
approximation to the Hashin model was utilized. However, in the future study of impact 
problems, where the out-of-plane stresses will be significant, three-dimensional versions of the 
criterion are available and will be used [50]. 
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The Hashin criteria were based on quadratic combinations of stresses and strengths. 
Quadratic equations were chosen in order to provide the best approximation to the failure surface 
while still allowing for a relatively simple model [50], Ply failure based on fiber tensile failure, 
fiber compressive failure, matrix tensile failure and matrix compressive failure was predicted 
separately. In each of the separate criteria, failure was considered to have occurred if the value of 
the expression is greater than one (1). For the purposes of this study, once failure in any of the 
failure modes was detected, total composite failure was considered to have occurred. In actuality, 
particularly for matrix dominated failure modes, the composite could withstand load after the 
“failure” load has occurred. However, for the simple uniaxial off-axis laminate configurations 
examined in this study, final composite failure was assumed to occur shortly after initial matrix 
cracking takes place, therefore only the initial failure load was predicted. 

Failure criteria for each of the local failure modes were as follows. In each of the 
expressions, Gjj was the macroscopic stress component, X T was the ply tensile strength in the 
longitudinal (fiber) direction, and X c was the compressive strength in the longitudinal direction. 
Furthermore, Y T was the tensile strength in the transverse direction (perpendicular to the fiber in 
the plane of the composite), Y c was the compressive strength in the transverse direction, and X s 
was the ply in-plane shear strength. Failure was assumed to occur when the value of the strength 
expression became greater than or equal to one (1). Tensile fiber failure was predicted by using 
the following expression: 



7 

'<7,2 > 

( X r i 

1 

k J 


= 1 


(5.1) 


Compressive fiber failure was predicted using the following equation. Shear stresses were not 
included in the failure criterion since Hashin was unsure whether shear stresses increased or 
decreased the compressive strength. Therefore, the effects of shear stresses were neglected [50]. 



Tensile matrix failure was predicted using the following expression: 
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5.3 Verification Analyses 


In order to verify the capabilities of the failure criteria as implemented into the composite 
micromechanics equations, the IM7/977-2 and AS4/PEEK composites considered in the last 
chapter were analyzed. The ply failure stresses were predicted for various fiber orientation angles 
and strain rates. 

5.3.1 IM7/977-2 Composite 

To verify the ply strength model, the IM7/977-2 composite considered in the last chapter was 
analyzed. For the IM7/977-2 system, the longitudinal tensile strength was 2300 MPa [67], the 
longitudinal compressive strength was 900 MPa [80], the transverse tensile strength was 73 MPa 
[67] and the shear strength was 85 MPa [67]. Due to a lack of data, the transverse compressive 
strength was assumed to be twice the transverse tensile strength. To compute the in-plane shear 
strength, the failure stress of a [±45°] s laminate was divided by two (2), which is a standard 
procedure for determining shear strength [78]. The predicted and experimental [67] failure 
strength values for [10°] and [45°] laminates for the IM7/977-2 material system are shown in 
Table 5.1. Note that the experimental values of the failure stresses likely have some scatter, but 
the statistical data were not available. For both laminates considered, failure was predicted to be 
due to tensile matrix failure. The “failure stress” stated in Table 5.1 and all remaining tables was 
the longitudinal stress in the structural axis system (along the loading direction) at which failure 
was predicted to occur. However, all of the stress components were used in applying the Hashin 
failure criteria. 

As can be seen in Table 5.1, for both fiber orientations the predicted failure stresses 
compared favorably to the experimental values, and most likely fell within the experimental 
scatter. These results indicated that the presented failure criteria produced accurate results for a 
variety of fiber orientations. However, more strength tests should be conducted over a wider 
range of strain rates to fully verify the model. 

5.3.2 AS4/PEEK Composite 

The ply strengths of the AS4/PEEK composite studied in the last chapter were also predicted. 
For this material, only quasi-static strength data were available. Ply shear strength data that 
provided an acceptable correlation with the available experimental results were not available. 
Furthermore, transverse stresses predicted using the deformation model for off-axis composite 
layers (such as [30°] and [45°] laminates) were greater than the transverse strengths indicated by 
the Weeks and Sun [25] data shown in Figures 4.9 and 4.13. Therefore, these values were not 
used for the transverse strengths. 

Figures 4.9 and 4.13 indicated that the transverse modulus did not appear to vary with strain 
rate for this material, indicating that the transverse strengths were also rate independent. For a 
carbon fiber reinforced polymer matrix composite, longitudinal strengths have been found to be 
rate independent [7]. Therefore, for this study, the longitudinal and transverse strengths for the 
AS4/PEEK system were assumed to be rate independent. However, as can be seen from Figures 
4.6-4. 8 and Figures 4.10-4.12, for off-axis fiber orientations the strengths appeared to be rate 
dependent. Therefore, for this study, the shear strength was assumed to be rate dependent and the 
cause of the rate dependence seen in the strengths of the off-axis laminates. 
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For the AS4/PEEK material, a longitudinal tensile strength of 2070 MPa was used [81], and 
the transverse tensile strength was set to 83 MPa [81]. The longitudinal compressive strength 
was set equal to one-half of the longitudinal tensile strength, and the transverse compressive 
strength was again assumed equal to twice the transverse tensile strength, based on representative 
composite data [78], The shear strength values were computed by correlating results from the 
[15°] laminates (Figures 4.6 and 4.8). From this information, the shear strength for a strain rate 
of lxlO' 5 /sec was determined to be 63 MPa, and the shear strength for a strain rate of 0.1 /sec 
was determined to be 88 MPa. All of these strength values likely had some scatter, but the 
statistical data were not available. 

Failure stresses were predicted for the [30°] and [45°] laminates for strain rates of lxlO' 5 /sec 
and 0.1 /sec. The predicted and experimental results for a strain rate of lxlO' 5 /sec are shown in 
Table 5.2, and the results for a strain rate of 0.1 /sec are shown in Table 5.3. In all cases, failure 
was predicted to be due to tensile matrix failure. 

For both strain rates and both fiber orientations considered, the comparison between the 
predicted and experimental values was quite good and most likely within the experimental 
scatter. The results indicated that the failure criteria were able to predict ply failure for a variety 
of fiber orientations and strain rates. The results for AS4/PEEK also indicated that even when 
some approximations were required in determining the ply failure stresses, reasonable results 
could still be obtained. 

5.4 Summary 

The Hashin failure criteria were implemented into the composite micromechanics equations 
in order to predict the ply failure stresses in polymer matrix composites. The failure criteria 
predicted ply ultimate strengths by using quadratic combinations of macroscopic stresses and ply 
strengths to approximate local failure mechanisms. Equations were available to predict ply 
failure based on fiber tensile or compressive failure, or matrix tensile or compressive failure. 
Verification studies were conducted using two representative polymer matrix composites. The 
results for the IM7/977-2 system showed that if good macroscopic strength data can be obtained, 
ply ultimate strengths could be reasonably predicted with no approximations required. The 
results for the AS4/PEEK material, on the other hand, demonstrated that even if some 
approximations in the ply failure strengths were required, reasonable results could still be 
obtained. The AS4/PEEK calculations also provided preliminary indications that the rate 
dependence of ply strengths can be predicted. However, final application of the Hashin model 
will require developing expressions for determining the transverse and shear strengths as a 
function of strain rate. Strength tests conducted over a wide range of strain rates will help to 
accomplish this goal. The predictions for both material systems showed that the ply ultimate 
strength could be accurately determined for a variety of fiber orientation angles. 
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Table 5.1 

Failure Stress Predictions for IM7/977-2 Laminate 



Predicted Failure Stress 
(MPa) 

Experimental Failure Stress 
(MPa) 

[10°] Laminate 

480 

500 

[45°] Laminate 

100 

105 


Table 5.2 

Failure Stress Predictions for AS4/PEEK at Strain Rate of lx 10' 5 /sec 



Predicted Failure Stress 
(MPa) 

Experimental Failure Stress 
(MPa) 

[30°] Laminate 

130 

140 

[45°] Laminate 

98 

104 


Table 5.3 

Failure Stress Predictions for AS4/PEEK at Strain Rate of 0.1 /sec 



Predicted Failure Stress 
(MPa) 

Experimental Failure Stress 
(MPa) 

[30°] Laminate 

165 

170 

[45°] Laminate 

114 

112 
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CHAPTER 6 

FINITE ELEMENT IMPLEMENTATION 


The polymer constitutive equations and the composite micromechanics model described in 
earlier chapters were implemented into LS-DYNA, a commercially available transient dynamic 
finite element code. In this chapter, an overview of LS-DYNA and the user defined material 
option will be given. The transformation of the polymer constitutive equations and the composite 
micromechanics model into an incremental form will be described. Verification studies 
conducted using the PEEK thermoplastic and the AS4/PEEK composite discussed in previous 
chapters will be presented. 

6.1 Overview 


In order to simulate the impact response of a composite structure, the ultimate long term goal 
of this research, finite element methods are required. With that goal in mind, the matrix 
constitutive equations and the composite micromechanics model described in previous chapters 
were implemented into LS-DYNA, a commercially available transient dynamic finite element 
code [82]. LS-DYNA uses explicit central difference integration methods to integrate the rate 
equations. 

LS-DYNA had several material options for the analysis of composite materials. However, all 
of these models were for strain rate independent material response. Furthermore, only limited 
abilities to model nonlinearities in shear were present in the models, based on the Chang and 
Chang model [52,53]. 

LS-DYNA had a user defined material option that allowed users to implement material 
models that were not included with the finite element package. To use this option, a FORTRAN 
subroutine would be created which could then be linked to the finite element code. In the 
subroutine, strain increments for the current time step were passed in from the main program. 
Stresses and other history variables (such as state variables) from previous time steps were also 
available to the subroutine. From this information, the user defined material subroutine computed 
the stress increments, total stresses, and values of the history variables at the end of the current 
time increment. 

6.2 Incremental Form of Polymer Constitutive Equations 

To implement the polymer constitutive model discussed earlier into an LS-DYNA user 
defined material subroutine, the equations were converted into an incremental format. 
Specifically, Equations (3.7), (3.19) and (3.20) were adjusted so that increments in inelastic 
strain, internal stress and effective inelastic strain were computed instead of their corresponding 
rates. Rocca and Sherwood [62] implemented a version of the Ramaswamy-Stouffer constitutive 
equations into LS-DYNA to analyze the rate dependent inelastic deformation of polycarbonate. 
Their methodology was used as the basis for implementing the constitutive model described in 
this report. Rocca and Sherwood designed their user designed material subroutine to be used 
with plane stress shell elements. The same plane stress assumptions were used in the model 
implementation discussed here to maintain consistency with the implementation of the composite 
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micromechanics equations discussed later in this report. Maintaining the plane stress assumption 
for the polymer constitutive equations facilitated verification of the model. 

To convert the flow equation (Equation (3.7)) into an incremental form, the rate equation was 
multiplied by the time increment At of the current time step to compute the inelastic strain 
increment Ae,y . The resulting equation was as follows: 


Ae,' = D„ exp 
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where all of the terms were as defined earlier. Note that the total value of the deviatoric stress 
components and the internal stress components were used instead of the stress increments, and 
were the values from the previous time step. Equations (3.19) and (3.20) were modified to 
compute the increment in internal stress, AQjj, and the increment in effective inelastic strain, 
A£e'. The following equations resulted: 


aq,,=-^,,A£'-<A,A£; 
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and all the terms in the equations were as defined earlier. In Equation (6.2), the total value of the 
back stress from the previous time increment was used in computing the stress increment for the 
current time step. 

By examining Equations (6.1)-(6.3), one can notice that in many respects the transformation 
of the rate equations into an incremental format was carried out by applying forward Euler 
numerical integration methods [65]. While very simple equations resulted from this process, the 
forward Euler method would not be as accurate as other methods, particularly for stiff 
differential equations. As will be discussed in the analysis section of this chapter, using this 
simple approximation to integrate the rate equations lead to some discrepancies in the numerical 
results. 


6.3 Numerical Implementation of Polymer Constitutive Equations 

The incremental form of the polymer constitutive equations was incorporated into LS-DYNA 
through the use of a user defined material subroutine. As a first step, only the polymer 
constitutive equations were implemented and tested. The composite micromechanics equations 
were added once the constitutive model was verified to be programmed in correctly. 

In the user defined subroutine, the stresses, deviatoric stresses and internal stresses from the 
previous time step, along with the strain increment for the current time step, were passed into the 
routine. The material properties were also read in from the LS-DYNA input file. The inelastic 
strain components for the current time step were computed by taking the inelastic strain rate 
determined in the previous time step, and multiplying it by the current time increment, as 
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suggested by Equation (6.1). The effective inelastic strain increment was computed using 
Equation (6.3). 

For the shell element formulation, the through-the-thickness strain increment Ae.n was 
computed using the following expression: 

Ae„ = -j— — (ae/; + Ae 22 )+ Ac', (6.4) 

where v was the Poisson’s ratio and A£jj E was the elastic component of the strain increment. 
From this information, the mean and deviatoric strain increments were computed using standard 
definitions [24], The elastic component of the deviatoric strain increment was determined by 
assuming that the inelastic portion of the deviatoric strain increment was equal to the inelastic 
strain increment [62], This assumption neglected volume effects, which might be significant for 
the polymers considered in this study. Once the strain increments were computed, the stress 
increments and the deviatoric stress increments were calculated. The total stresses and total 
deviatoric stresses for the current time step were computed by adding the stress increments to the 
total stresses from the previous time step. 

The internal stress increment was computed using Equation (6.2), and this value was added 
to the total internal stress from the previous time step to determine the total internal stress for the 
current time step. Equations (3.7) and (3.9-3.18) were used to compute the inelastic strain rate 
components for the current time step. The inelastic strain rates, deviatoric stresses and internal 
stresses were then stored in history variables to be used in future time steps. 

6.4 Verification Analyses for Polymer Constitutive Equations 

To test the implementation of the polymer constitutive equations, the tensile response of the 
PEEK material analyzed in previous chapters was computed. The elastic and inelastic material 
properties for PEEK are given in Table 3.2. The tensile stress-strain curves at constant strain 
rates of 0.1 /sec and 1.0 /sec were simulated using LS-DYNA. While these strain rates were 
slightly above the strain rate range for which the material was characterized, they were 
considered to be low enough for the material constants to still be valid. However, LS-DYNA was 
designed to be utilized for high strain rate applications, and sometimes had stability and 
convergence problems at lower strain rates [82], The strain rates used for these analyses were 
chosen to be high enough to hopefully minimize these difficulties, but low enough to ensure 
validity of the material constants. The relatively low strain rates used for the analyses could still 
be a cause for any discrepancies in the computed results, as the integrator might not work well 
for low strain rates. 

For the finite element model, four noded shell elements were used in a square mesh, ten 
elements on a side as is shown in Figure 6.1. Each side of the model was 20 mm long. The left 
hand side of the model was clamped, and a constantly increasing specified displacement was 
applied to the right hand side of the model. These boundary conditions were chosen in order to 
simulate a constant strain rate tensile test. The displacements applied at each time were 
computed by taking the constant strain rate, multiplying it by the current time to obtain total 
strain, and multiplying this value by the total length of the model to compute an average 
displacement. The finite element model was designed to simulate the behavior of the polymer at 
an infinitesimal material point. A model of this type was used instead of attempting to model the 
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actual specimen in order to facilitate comparison of the finite element results to the results 
computed using stand alone computer codes. A mesh convergence study showed that the mesh 
was adequate. To generate a stress-strain curve, results were read in from several elements at the 
center of model, in order to avoid edge effects. 

Stress-strain curves computed at a strain rate of 0.1 /sec are shown in Figure 6.2, and tensile 
curves computed at a strain rate of 1.0 /sec are shown in Figure 6.3. In the figures, results 
computed using LS-DYNA (labeled “Finite Element” in the figures) were compared to results 
computed using the stand-alone computer code described in Chapter 3 (labeled “Analytical” in 
the figures). Note, however, that since the strain rates considered here were higher than those 
examined in the earlier analyses, the computed stresses were higher than those shown in Figures 
3.6-3. 8. Since all of the elements in the center of the finite element model had approximately 
equal stresses, only one curve is shown for the finite element results. 

As can be seen in the figures, for both strain rates the finite element results compared 
reasonably well to the results computed using the stand-alone computer code. The nonlinearity 
and rate dependence of the tensile curves was captured by the finite element analysis. The finite 
element computations predicted somewhat higher stresses than the stand-alone code, particularly 
in the inelastic range. There were two possible causes for this discrepancy. First of all, as 
mentioned earlier the conversion of the constitutive equations to an incremental form used 
forward Euler types of approximations, which were less accurate than other types of integration 
methods, particularly for the stiff equations used here. In addition, the strain rates at which the 
analyses were conducted were significantly lower than is usually recommended for an LS- 
DYNA analysis. These factors most likely combined to cause the discrepancies seen in the 
analyses. As will be discussed in a later section of this report, these effects were most likely even 
more significant when the polymer constitutive equations were implemented into the composite 
micromechanics. Further discussion of the causes of the discrepancies will be presented in the 
section describing the verification analyses for the LS-DYNA implementation of the composite 
micromechanics equations. 

6.5 Incremental Form of Composite Micromechanics Equations 

The general format of the micromechanics equations as implemented into LS-DYNA were 
the same as discussed in Chapter 4. However, one difference was that in the user defined 
material subroutine stress increments were computed in each subcell using strain increments and 
inelastic strain increments. This is different from the stand-alone code described earlier, in which 
total stresses were computed based on total strains and total inelastic strains. Furthermore, for the 
LS-DYNA implementation a plane stress condition was assumed on the macroscopic unit cell 
level, and shell elements were used in the finite element analyses. The plane stress condition was 
assumed since the use of shell elements in LS-DYNA simplified the analysis of laminated 
composite material as compared to using solid elements. However, within the subcells of the 
composite unit cell, a full three-dimensional stress state was still permitted. The macroscopic 
through-the-thickness strain increment, A£ 33 , was computed using the following expression [36]: 

Ae,_, = 11 (ae, , - Ac/, ) 21 (a£, 2 - Ae ' 22 )+ Ac,, (6.5) 

C 3? Cjj 
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where A£jj was the strain increment, Ae^ 1 was the effective inelastic strain increment, and Qj 
represented terms in the macroscopic elastic stiffness matrix. The relationship between the 
elastic stiffness matrix components and engineering constants such as modulus and Poisson's 
ratio was found in standard texts such as [36.47]. The effective inelastic strain increment was 
determined by multiplying the effective inelastic strain rate by the value of the time increment 
At, similar to what was described in Equation 6. 1 . 

6.6 Numerical Implementation of Composite Micromechanics Equations 

The incremental form of the composite micromechanics equations was once again 
implemented into LS-DYNA through the use of a user defined material subroutine. In the user 
defined routine, the macroscopic strain increments for the current time step, and the stresses and 
internal stresses for each subcell from the previous time step, were passed into the routine, along 
with the material properties. The through-the-thickness strain increment was computed using 
Equation 6.5. The inelastic strain increments were computed by multiplying the inelastic strain 
rates by the value of the time increment. 

The stress increments and total stresses for each subcell were computed using the composite 
micromechanics equations, and the total stresses for the unit cell were determined using the 
uniform stress and uniform strain assumptions described earlier. The inelastic strain rate and 
internal stress for each subcell was determined using the incremental form of the polymer 
constitutive equations. The effective inelastic strain rates for the unit cell were calculated using 
Equations (4.59)-(4.64). The inelastic strain rates, internal stresses and total stresses for each 
subcell were stored in history variables for use in future time steps. 

6.7 Verification Analyses for Composite Micromechanics Equations 

The implementation of the composite micromechanics equations was tested by analyzing the 
AS4/PEEK composite considered in Chapter 4. Unidirectional composites with fiber orientations 
of [15°], [30°], [45°] and [90°] were once again considered. A constant strain rate of 0.1 /sec was 
applied to the finite element model. This strain rate was once again significantly lower than was 
usually recommended for an LS-DYNA analysis. However, this strain rate was within the range 
for which the PEEK material was characterized, and there were experimental data available at 
this strain rate. The material properties of the AS-4 fibers are given in Table 4.1, and the elastic 
and inelastic properties of the PEEK thermoplastic are once again given in Table 3.2. 

The finite element model shown in Figure 6. 1 was once again used. The boundary conditions 
described in the analysis of the pure PEEK material were applied to the model. An important 
point to note is that for the off-axis fiber orientation angles, the boundary conditions did not truly 
represent the periodicity conditions required to simulate the deformation response of a material 
point. However, periodic boundary conditions would have been overly difficult to apply to the 
model. The fiber orientation angle was specified in the shell element definition of the LS-DYNA 
input deck. To generate the stress-strain curves, results were read in from elements 45-47 and 55- 
57 of the model (see Figure 6.1). These elements, near the center of the model, were chosen in 
order to minimize edge effects and to minimize to the greatest extent possible the effects of the 
boundary conditions. 
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The computed stress-strain curves are shown in Figures 6.4-6. 7. In each figure, results 
predicted by the finite element analyses were compared to the experimental values obtained by 
Weeks and Sun [25]. For the [30°] and [45°] laminates, the predicted finite element results varied 
significantly from element to element, even in the center of the model. This variation was most 
likely due to the nature of the approximate boundary conditions that were applied. Therefore, 
results from several elements were plotted in the figures. The elements that were used to generate 
each curve were listed in the legend of each figure. The results from the elements that are not 
plotted fell within the bounds given by the plotted curves. The actual tensile response, if proper 
boundary conditions were applied, would most likely have fallen within the bounds of the plotted 
curves. 

As can be seen in the figures, for all of the fiber orientations the elastic portion of the 
deformation response was predicted quite well by the finite element analyses. However, the 
stresses in the inelastic range were overpredicted when compared to the experimental values. 
Furthermore, as can be seen in the results for the [15°] and [30°] laminates, once saturation was 
about to occur, the finite element results became somewhat unstable, with the total stresses 
beginning to rise with increasing strain. 

There were several possible causes for the discrepancies seen in the predictions of the 
inelastic stresses. First of all, the fact that the applied boundary conditions did not truly represent 
proper periodic boundary condition most likely affected the results. For the off-axis fiber 
orientations, the response of the finite element model was not truly representative of the response 
of an infinitesimal material point. In addition, since the value of the material constant (3 was 
determined empirically, some modification of this value might be required. 

More significant explanations for the discrepancies between the experimental and predicted 
results, however, related to the nature of the incremental formulation of the polymer constitutive 
equations. As discussed earlier, in the analyses of the pure PEEK thermoplastic the inelastic 
stresses were somewhat overpredicted. The reasons for those discrepancies were hypothesized to 
be due to both the low strain rate which was used for the analysis and the method used to 
integrate the rate equations in the polymer constitutive model. Furthermore, in the analyses of the 
pure PEEK material, an essentially uniaxial stress state was applied to the material. For the 
composite, on the other hand, the polymer was in a fully multiaxial stress state, which might 
have magnified any discrepancies and instabilities in the analysis. 

To further explore the effects of the integration method used on the rate equations in the 
polymer constitutive model, the implementation of the equations was revised. Specifically, each 
time step was broken down into ten equally spaced “substeps” within the user defined material 
subroutine. The Euler integration was then carried out for each “substep” in the analysis. By 
reducing the size of the time step over which each integration was performed, the accuracy of the 
computations should improve. Specifically, the methodology described in Section 6.6 was again 
utilized. However, each calculation was carried out ten times for each time step, using time and 
strain increments equal to one-tenth of the time and strain increments passed into the routine. 

Finite element analyses were carried out on the [15°] and [30°] AS4/PEEK laminates at a 
strain rate of 0. 1 /sec using the revised equation formulation, and the results are shown in Figures 
6.8 and 6.9. By comparing Figure 6.8 to Figure 6.4, and by comparing Figure 6.9 to Figure 6.5, 
the observation can be made that revising the method used to integrate the polymer constitutive 
equations improved the results both qualitatively and quantitatively. These results indicated that 
the method used to integrate the rate equations in the finite element implementation could 
significantly affect the results. However, since the revised algorithm significantly increased the 
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execution time of the analysis, more sophisticated integration methods need to be implemented 
in the future. However, by further refining the integration methods and using higher strain rates 
in the analyses, the discrepancies between the predicted and experimental results should be 
significantly reduced. 

6.8 Summary 


The polymer constitutive equations and the composite micromechanics model have been 
implemented into the explicit finite element code LS-DYNA. The analytical equations described 
in previous chapters were converted into an incremental format. Verification studies were 
conducted on the PEEK thermoplastic and the AS4/PEEK composite studied in previous 
chapters. The stresses predicted in the inelastic range of the deformation response were 
somewhat overpredicted compared to experimental results and results computed using stand- 
alone computer codes. Possible explanations for the discrepancies were identified and discussed. 
The method used to integrate the rate equations was found to have a significant effect on the 
accuracy of the results. Future work will concentrate on improving the integration methods used 
in converting the rate equations into an incremental form. In addition, methods of improving the 
boundary conditions applied when studying off-axis composite laminates will be considered. 
Furthermore, once high strain rate data are obtained, characterization studies and analyses will be 
conducted using this data. Since LS-DYNA was designed to be used at higher strain rates, 
analyses conducted at high strain should compare more favorably to the experimental data. 
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Figure 6.1: Finite Element Model Used for Verification Analyses 
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Figure 6.2: Finite Element Predictions for PEEK at Strain Rate of 0.1 /sec 



Figure 6.3: Finite Element Predictions for PEEK at Strain Rate of 1.0 /sec 
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Figure 6.4: Finite Element Predictions for AS4/PEEK [15°] Laminate 

at Strain Rate of 0. 1 /sec 



Figure 6.5: Finite Element Predictions for AS4/PEEK [30°] Laminate 
at Strain Rate of 0. 1 /sec 
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Figure 6.8: Finite Element Predictions for AS4/PEEK [15°] Laminate 
at Strain Rate of 0. 1 /sec with Revised Equation Formulation 
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Figure 6.9: Finite Element Predictions for AS4/PEEK [30°] Laminate 
at Strain Rate of 0. 1 /sec with Revised Equation Formulation 
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CHAPTER 7 

SUMMARY AND CONCLUSIONS 


A rate dependent deformation and strength model for polymer matrix composites has been 
developed which accounts for material nonlinearities. The deformation model has also been 
implemented into a transient dynamic finite element code. The goal of this research was to 
develop preliminary models that would provide insight into the analytical methods required to 
conduct detailed simulations of polymer matrix composites subject to high strain rate impact. In 
this chapter, a general summary of the work presented in this report will be given, along with a 
discussion of future work. 

7.1 Summary 

Rate dependent, inelastic constitutive equations based on the Ramaswamy-Stouffer state 
variable equations were formulated and implemented numerically to model the nonlinear 
deformation of ductile, crystalline (or semi-crystalline) polymers. While the Ramaswamy- 
Stouffer equations were originally developed to simulate the response of metals above about 
one-half of the melting temperature, appropriate modifications were made to the model in order 
to analyze the response of polymers. For example, the effects of hydrostatic pressure on the 
inelastic response were accounted for by modifying the effective stress definition in the flow 
law. Specifically, the shear terms in the effective stress were modified to model the effect of 
hydrostatic stresses on the deformation response. Correlation studies were conducted using two 
representative polymeric materials, Fiberite 977-2 and PEEK, and good correlation between the 
computed deformation response and experimental values was obtained. 

A mechanics of materials based micromechanics technique was developed and numerically 
implemented to predict the effective deformation response of a composite with a nonlinear, rate 
dependent matrix constituent. The effective response of the composite was computed based on 
the response of the individual constituents. The fibers were assumed to be linearly elastic, and 
the polymer matrix was analyzed using the constitutive equations described above. To formulate 
the micromechanics equations, uniform stress and uniform strain assumptions were applied to a 
composite unit cell. The effective stresses in the composite unit cell were computed given the 
effective strains. Constituent level stresses in the fiber and matrix were also calculated, along 
with the effective inelastic strains. Verification studies were conducted using two polymer matrix 
composites, IM7/977-2 and AS4/PEEK. Analyses were conducted on these materials for a 
variety of fiber orientations and strain rates. The predicted results compared well to 
experimentally obtained stress-strain curves. 

To predict the ultimate strength of a composite ply, which is a necessary step in the 
development of full failure and penetration models for a polymer matrix composite under high 
strain rate impact, the Hashin failure criteria were implemented into the composite 
micromechanics equations. In the Hashin criteria, local failure mechanisms such as fiber failure 
or matrix cracking were approximated based on ply level stresses and strengths. Ply ultimate 
strengths were predicted for the two composites considered in the earlier analyses, and the results 
were compared to experiment. Even though some approximations were required in determining 
the ply strengths for the AS4/PEEK material, for both materials the comparison between the 
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predicted ply strengths and the experimental values were very good for a variety of fiber 
orientations and strain rates. 

The polymer constitutive equations and the composite micromechanics deformation model 
were implemented into LS-DYNA, a transient dynamic explicit finite element code. The user 
defined material subroutine option was used in implementing the equations into LS-DYNA. To 
conform to the requirements of LS-DYNA, the equations were converted into an incremental 
format. Given strain increments, increments in stress, internal stress and inelastic strain were 
computed. The implementation of the deformation model was carried out by analyzing the PEEK 
thermoplastic and AS4/PEEK composite considered earlier. The results computed for the pure 
polymer by finite element analysis compared reasonably well to stress-strain curves computed 
using a stand-alone computer code, with some overprediction of the stresses in the inelastic 
range of the deformation response. The stress-strain curves generated by finite element analysis 
for the AS4/PEEK composite compared quite well to experimental values in the elastic range. In 
the inelastic range, the stresses were generally overpredicted, and the analysis displayed some 
instability after “saturation” was reached. Further analyses indicated that the integration method 
used to implement the incremental form of the polymer constitutive equations likely was a 
significant cause of the discrepancies, as refining the integration technique led to improvements 
in the calculation of the inelastic stresses. The low strain rates at which the analyses were carried 
out most likely also contributed to inaccuracies in the finite element calculations. 

7.2 Future Work 


As mentioned earlier, the overall goal of this research was to develop preliminary 
deformation and strength models that would provide insight into the analytical methods required 
to conduct detailed simulations of polymer matrix composites subject to high strain rate impact. 
As a result of this study, several areas of research were indicated as future steps in accomplishing 
the overall goal. First of all, the polymer constitutive equations will be revised in order to 
represent the physical mechanisms more accurately, particularly the dependence of the 
deformation response on the hydrostatic stress. Furthermore, the equations will be modified so 
that the deformation during unloading can be predicted by the model. The constitutive model 
will also be extended into the large deformation domain. Finally, more sophisticated integration 
methods will be utilized in the numerical implementation of the model in order to improve the 
robustness of the calculations. 

A lamination theory will be incorporated into the composite micromechanics equations, to 
allow for the analysis of full composite laminates. The model will also be modified to allow for 
the analysis of woven composites, as these types of composites are more likely to be used in fan 
containment applications. Finally, the micromechanics will also be extended into the large 
deformation regime. 

The ply strength model will be extended to allow for the full penetration and failure analysis 
of composite laminates. A property degradation model will also be incorporated into the failure 
calculations in order to allow for the gradual reduction of material properties during impact 
penetration. In addition, the developed methods will also be modified to allow for the analysis of 
woven composites. 
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All of the modifications mentioned above will be implemented into LS-DYNA through 
modification of the user defined material subroutine. Furthermore, the methods used to develop 
incremental forms of the constitutive equations will be modified in order to improve the stability 
and accuracy of the finite element calculations. The developed penetration and failure models 
will also be implemented into LS-DYNA as part of this process. 

High strain rate Hopkinson bar tests will be conducted on representative materials, and the 
constitutive equations will be characterized and validated for high strain rate conditions. High 
strain rate impact tests will then be conducted on laminated and woven polymer matrix 
composites. The impact tests will then be simulated using finite element analyses. The ultimate 
goal of this research is to be able to accurately simulate the high strain rate impact of polymer 
matrix composites. Once this goal is accomplished, new and improved design tools will thus be 
available which can help engineers design structures subject to high rate loading such as aircraft 
engine fan containment systems. 
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